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MODIFICATIONS OF THE CYGRO-2 COMPUTER PROGRAM 

by Ivan B. Fiero 
Lewis Research Center 

SUMMARY 

The CYGRO-2 computer program is a potentially useful tool in analytically predict- 
ing fuel-element swelling for the cylindrical geometry, nuclear reactor fuel pins. The 
program takes into account many facets of the behavior of irradiated fuel pins. These 
include fuel swelling, fuel-clad interaction,, thermal gradients, and plastic and creep 
behavior. However, the program as it existed was not compatible with the Lewis 
IBM 7094 computer. It had inherent problems associated with application to this com- 
puter and associated with techniques used to solve equations within the program. 

The program was modified such that it would be compatible. However, many diffi- 
culties were encountered while using the program. Termination before completion often 
occurred in conjunction with various sorts of computer diagnostics. These computer 
diagnostics, once loading was achieved, consisted mainly of overflows, underflows, 
divide by zero, and too large exponents. Also, the program often utilized very small 
time steps . Several modifications were made to the program to alleviate these diffi- 
culties. These modifications can be summarized as three main categories: (1) modi- 
fications in the techniques of solving certain simultaneous equations and methods of 
averaging variables for a given time step, (2) optional changes such as printout options 
to aid in debugging and some additional input and output which are useful to the program 
user, and (3) changes in FORTRAN form and range of subscripts and general shortening 
for application to the Lewis computer. 

The resulting modifications not only caused the program to run as was primarily 
desired, but resulted in decreased computation times and improved convergence. The 
improved computation times were very significant in some instances. 

A sample problem of a uranium mononitride fuel pin has been included to show the 
changes that have been made to the input data deck and the printout received from the 
modified CYGRO-2. Changes to the data deck include modification of the initial com- 
ments cards and an inclusion of two additional pieces of data. The printout format was 
changed slightly and a few more variables were added to the print list. Various prob- 
lems associated with the input data are also discussed. 



INTRODUCTION 


The CYGRO-2 computer program (ref. 1) is a potentially useful tool in analytically 
predicting fuel-element swelling and growths for the cylindrical geometry, nuclear 
reactor fuel pins. It was written by Friedrich and Guilinger of Westinghouse- Bettis as 
a development from a group of individual computer programs (such as FIGRO, CYGRO, 
and CYGRO- 1; refs. 2 to 6) to calculate the effect of swelling fuel, fuel- clad interac- 
tions, and pin design. Basically, the program is based on a two-material zone (fuel 
zone and clad zone), cylindrical geometry fuel pin constrained by a plane strain assump- 
tion. This assumption best applies to long cylinders, where the effect of end constraints 
and multiple pellets are neglected. 

The fuel may or may not contain a central void or plug and there may be an initial 
gap between fuel and clad. The program utilizes viscoelastic behavior determined by 
both a normal creep rate associated with second- stage creep and an accelerated creep 
which corresponds to plastic flow. Usual constraints of equilibrium and compatibility 
are assumed. Mechanical loads originate from thermal expansions, gas pressures, and 
irradiation swelling. Irradiation swelling depends on (1) solid fission products gener- 
ated and deposited within the fuel and fragmented into the clad and (2) gaseous fission 
products. All the swelling associated with gaseous fission products is assumed to be 
the result of coalescing into a set number of inter- and intragranular bubbles. Gas 
bubble growth depends on an equilibrium between bubble pressure, surface tension, 
creep forces adjacent to the bubble, and the spherical stress in the surrounding mater- 
ial. Stress in the surrounding material is dependent on the forces arising from thermal 
gradients and fuel- clad interactions. Fuel cracking and fuel- clad interface slip are not 
considered in this version of CYGRO, although they are included in a later version 
(ref. 7). 

The program achieves solution by a finite difference formulation. Due to axial 
symmetry, the only geometric variable is radius, which is divided into thin cylinders 
or rings. Also, the program is quasi-transient, which means small time increments 
or steps are utilized, during which rate variables are assumed to be constant. 

The advantage of such a complicated program as CYGRO-2 is that it takes into 
account much of the behavior that is known to take place within irradiated fuel pins. 

For example, it includes thermal expansions, creep and plastic flow, generation of 
fission products requiring more volume than the initial material occupied, releasing of 
gas into plenums or voids, and sintering and hot pressing behavior. However, there are 
certain disadvantages in that (1) input is lengthy and complex, (2) some input is difficult 
to evaluate (such as gas gap pressure, which can be strongly dependent on fuel swelling), 
(3) the complex interrelations of parameters involved makes it difficult to predict the 
effect of changes in certain parameters in advance, and (4) such a complex program 
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requires considerable computer time. A direct disadvantage of the complex input is 
that, in order to generate proper data, knowledge of fuel behavior is required, a priori. 
For example, the solid volume expansion of bare fuel must be known and the number of 
bubbles as a function of temperature must be known. These are only roughly or empir- 
ically known from previous fuel-pin tests, or are only known for some particular final 
burnup at the present stage of understanding of swelling of nuclear fuels. Hence, even 
though the program is versatile it is left to the discretion of the user to model the fuel 
pin adequately. 

Several difficulties were encountered with the use of CYGRO-2 on the Lewis 
IBM 7094 computer. These consisted at the outset mostly of underflows, overflows, 
divide by zero, too large exponents, and basic core storage problems. In addition, the 
program was utilizing very small time steps, resulting in large running times for small 
problem times. These operating difficulties occurred from time to time and were de- 
pendent on and sensitive to fuel-pin geometry and power history. 

Resulting modifications made to the CYGRO-2 computer program were in three 
main categories: (1) modifications in the techniques of solving certain simultaneous 
equations and methods of averaging variables for a given time step; (2) optional changes, 
such as printout options, to aid in debugging and some additional input and output; and 
(3) initial changes in FORTRAN form and range of subscripts and, finally, general 
shortening for application to the Lewis computer. These modifications not only caused 
the program to run as was primarily desired, but also resulted in decreased computa- 
tion times. The improved computation times were very significant in some instances. 
Another advantage of some of the modifications was an improved accuracy. 

This report summarizes the modifications made to CYGRO-2 and presents a sample 
problem which was solved using the modified version of CYGRO-2. No attempt is made 
to discuss the modifications on a FORTRAN statement basis, but rather in general 
terms, with an attempt to be as consistent as possible with symbols and terminology of 
reference 1. It is assumed herein that the reader is familiar with reference 1. 

The modifications made to the program are presented first. Then a detailed sample 
problem for a uranium mononitride fuel pin follows, with a discussion of the input data 
cards. Included under the data card description is a discussion of difficulties with the 
input data in relation to results obtained on the computer. 
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MODIFICATIONS 


Calculation Technique 

Conductivity as a function of porosity . - The program requires conductivity for the 
fuel and clad in order to determine temperature distributions within the fuel pin. The 
conductivity can be expressed as a function of the porosity or density of the material by 
appropriate constants, as discussed in reference 1. 

An error can occur, however, if only a fraction of the solid material, rather than 
100 percent, is associated with the bubbles. If this is done, the apparent density within 
the program is not the true density and the resulting conductivity is incorrect. This is 
due to the particular parameter being utilized to make the density corrections. In fact, 
the solid volume fraction unaccounted for, or not utilized around any type of bubble, 
becomes part of the porosity. The equations were modified slightly so that the gas 
bubble volume is summed over all bubble types. The sum is divided by the total gas- 
plus-solid volume to obtain a density ratio V^/V. This is the ratio described in ref- 
erence 1 . One of the FORTRAN common blocks had to be moved into the subroutine to 
make the correct variables common for the modification. These density corrections are 
performed in subroutine SEGT. 

Variable gas generation rate in fuel rings . - Required gas bubble information for 
the fuel pin includes the moles of initial gas, fission gas generation rate coefficients, 
diffusion or leakage coefficients, volumes of the bubbles and solid material volumes 
associated with the bubbles for each type of bubble, and the number of bubbles. The 
only variables dependent on both the material, fuel or clad, and the individual rings are 
the moles of initial gas, the volumes of gas bubbles, and the number of bubbles. The 
others are only material dependent. However, this specified information is for each 
bubble type and the dependency on rings and/or material can be specified for each bubble 
type. 

In some instances, it was found that it would be advantageous to have the fission gas 
generation coefficient as a function of the rings. Consequently, the subscript was" 
changed from a dual subscript on the gas bubble types and material to a subscript set on 
the gas bubble types and rings Both subscript variable FORTRAN name and subscript 
range dimension required changing. In this way a radially dependent gas generation can 
be utilized, if desired, which would reflect a radial flux variation. These changes were 
made in the subroutine dealing with gas bubble calculations, subroutine SGAS. 

Gas bubble equilibrium . - Bubble equilibrium and bubble expansion or contraction 
are determined in the subroutine SGAS. The gas bubble is acted upon by several forces: 
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(1) Gas pressure due to initial gas present (fabrication porosity) and to fission 
gases which accumulate as bubbles, (Buildup of gas pressure tends to expand the bubble 
and cause the fuel growth phenomenon. ) 

(2) Surface tension of the gas bubble - solid interface, which tends to reduce or re- 
strict bubble growth. (This is the sintering force (ref. 1) which, under application of no 
other force, tends to sinter the material from its fabricated porosity. ) 

(3) Macroscopic stress in the surrounding material, which acts as a hydrostatic 
pressure and can cause the hot pressing of partially dense material (ref. 1). 

(4) The equivalent creep stress set up in the material due to a change in bubble size 
(creep rate is directly related to the bubble growth rate), which also acts to restrict 
bubble growth. 

A summation of these forces must be zero at the bubble-solid interface for the bubble to 
be in equilibrium. 

The CYGRO-2 program, it was discovered, spent a major portion of the computer 
time doing iterations to converge on an equilibrium (balance of forces) situation. The 
technique of convergence consists of summation of forces against bubble growth rate 
until the summation is zero. It was not unusual for the program to take 400 to 500 iter- 
ations through a long series of equations to converge on a solution. This part of the 
program was rewritten to change the method of convergence. A typical case now con- 
verges in about six to eight iterations. This modification reduced overall run time by a 
factor of 4 or 5 for some problems. 

Previously, if an iteration resulted in nonconvergence in equilibrium, the next iter- 
ation merely adjusted the bubble growth rate by some very small amount and recycled 
through the calculations. Now, the new bubble growth rate, after two iterations over a 
much wider range of bubble growth rate, is modified according to slope and/or sign of 
the resultant force (or pressure) around the bubble (modified method of succession 
approximations, ref. 8). This is a rather complicated set of equations because the 
curves (bubble growth rate against summation of pressure forces at bubble-solid inter- 
face) describing the iteration scheme are not well behaved. In general, they are some- 
what hyperbolic and can have blips or small changes in slope along the asymptotic por- 
tions. In addition, three modes of creep or plastic flow can cause discontinuities. 

These three modes are associated with the physical problem in that the material sur- 
rounding the bubble may be only in creep flow, in both creep and plastic flow, or entire- 
ly in plastic flow. 

It was also found that occasionally the iteration scheme could not find an equilibrium 
condition for the bubble . An example of this is when the hyperbolic-shaped curve does 
not have a solution but becomes asymptotic to the axis representing a solution, or con- 
tains a small blip so that the curve appears to diverge from any solution. Previously, 
when this occurred, the answer closest to the equilibrium condition was used. It was 
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found that utilizing this answer results in an instability in the stress calculation for a 
given fuel node and also results in abnormal fuel swelling for that fuel node. These 
results, of course, affect the calculations in the time steps following. This particular 
problem would also exist with the new accelerated convergence method unless the pro- 
gram were modified such that the time step could be decreased. It was found that a 
smaller time step would always cause a shift in the curve, resulting in a solution. 

Consequently, the subroutine was modified with a routine that decreased the time 
step to a value such that bubble equilibrium could be achieved. The controlling variable 
in the routine is called JERK. Essentially, what happens is that when such a condition 
of nonequilibrium occurs, this variable (JERK) is set equal to the larger of two values, 
either 2 or a ratio proportional to how far from equilibrium the problem is. The pro- 
gram is then forced to exit from the subroutine SGAS for the particular node for which 
it was calculating gas bubble equilibrium. All the remaining nodes are cycled through 
to obtain the largest value of JERK. Then the program returns to subroutine STIM 
where the time step is divided by JERK and the calculations begin again. The number 
of times this nonequilibrium occurs is printed out at the end of the problem and has been 
as high as 25 in some of the problems run thus far. This is out of a thousand or more 
time steps. 

The routine was also later remodified to decrease the time step when two other 
problems arise. On occasion, the method of convergence set up would not be able to 
find a solution in a reasonable number of iterations. One occasion was when a small 
blip would occur. It is nearly impossible to detect these blips in all possible cases; and 
after several debugs to correct this situation, it was decided that a unique set of circum- 
stances could, on any new problem, result in not converging on a solution. Consequent- 
ly, after an unreasonable number of iterations (100 was deemed unreasonable) the pro- 
gram utilizes the JERK variable to reduce the time step. 

Secondly, there are three possible modes of creep and plastic behavior around the 
bubble, each of which greatly affects the shape of the curve. These are physically 
related to conditions where the material surrounding the bubble is in creep only, the 
material is partially in creep and partially in plastic flow, or the material is entirely in 
plastic flow. The program would often cycle between modes, looking for the correct 
solution. Occasionally, the cycling was excessive because the modes exhibited a dis- 
continuous behavior between modes; that is, one mode approached +°° , whereas the next 
mode approached -°° . Whenever the program cycled between modes an excessive num- 
ber of times, the time step was also decreased by the variable JERK. In this case, 
since the usual number of times was only once or twice, the program was set to utilize 
JERK if it cycled at least 25 times. Satisfactory results have been obtained for both of 
these situations, with no significant increase in the number of times JERK has been 
utilized. Apparently, this problem occurs only rarely. 


8 



These modifications to the subroutine SGAS have had the beneficial effect of both 
greatly decreasing computer time and improving accuracy. 

Fuel- clad gap closure . - Contact between the fuel and clad occurs because the radi- 
al growth of the fuel overtakes that of the clad. In the program the growth during a 
given time step may exceed that of the clad and give a sudden interference fit, providing 
the interference force does not exceed a specified amount. This interference is cor- 
rected in the next time step. However, the initial interference is frequently so great 
that it causes irregularity in the stresses and strains of the fuel and clad. Rates of 
stress and strain can be affected enough that, once contact occurs in this manner, the 
program calculates negative contact forces in the following time step. This period dur- 
ing the history of the fuel pin causes special problems and in the CYGRO-2 program 
(ref. 1) is called the ’’twilight" period. 

Now, rather than utilize the full time step, which would cause an interference fit/ 
the time step is proportionally reduced to cause gradual contact. A check is made to 
see if an interference is going to occur. If so, the time step is reduced an amount pro- 
portional to the ratio of the actual gap to the actual gap plus interference. An additional 
arbitrary multiplier of 0. 9 is used so that the contact is approached smoothly. This 
results in a few more iterations initially, but because of the smoother action resulting 
from this method of contact, the overall number of iterations is generally reduced. And 
there are no irregularities in the stresses and strains during this ’’twilight" period. 

Smallest time increment and gap closure criterion. - The CYGRO-2 computer pro- 
gram utilizes a very small number (10 ) to represent limitations on two parameters. 

It is defined by an arithmetic statement in the program and is alterable only by changing 
the FORTRAN. It is used (1) for the convergence criterion for the fuel- clad contact (the 
point at which the radii are set equal), and (2) as the smallest possible time step allowed 
(defined as TINY). The time step can decrease due to a number of factors: rates of 
stress or strain may be too large; fuel- clad contact interference may reduce the time 
step; or problems may arise with fission gas bubble equilibrium calculations, as pre- 
viously discussed. 

-20 

It was felt that 10 was too small for these purposes for two reasons. First, as 

a convergence criterion for the fuel- clad contact, it would be significant enough to use 

something much larger, say 10”^. Secondly, this small time step was not utilized by 

the program unless there were other problems which indirectly required a small time- 

step and needed corrective action separately. When the small time step (TINY) was 

used, none of the variables appeared to change. For example, it made no difference in 

-7 -20 

the values of the variables whether 10 was used or 10 . It simply took slightly 

longer to work down to the smaller number. Consequently, it was felt that the smaller 

-7 

number served no real useful purpose and it was changed to 10 . This means that, 

-4 -7 

when the gap is closed to within 10 mils (2. 54x10 cm), the radii are set equal. 


7 



These changes are required to the subroutine STIM and also in the main program to a re- 
lated variable called STINY. 

Another problem which has occurred, and is complicated by the previously dis- 
cussed modifications to subroutine SGAS, is that the program automatically bypasses 
the gas equilibrium calculations when this small time step (TINY) is used. It may be 
wise to further modify the program such that it cycles through the gas calculations even 
when the small time step is used. This is especially true during sintering when the 
bubble collapse results in the gas calculations being the governing factor for decreasing 
the time step. It is obvious that when the subroutine SGAS is responsible for the de- 
creased time step, it should not in turn be bypassed. This bypassing of subroutine SGAS 
is also the reason that sometimes none of the variables appear to be changing when 
TINY is used. A slight modification to the subroutine SGAS would cause the calculations 
to always be made. However, this change has not as yet been incorporated into the 
modified CYGRO-2. 

Smaller initial time step . - The initial time step, after a new time card was read, 
was originally chosen as 1/100 of the time interval between two consecutive time cards. 
If the problem required a lesser time, and it usually did, it could obtain it in several 
iterations by the routine previously discussed. This was modified to use the smaller of 
either 1/100 of the time interval or 0. 1 hour for the initial time step. This helps to 
reduce the number of iterations when there is a large interval of time between two con- 
secutive time cards. It does not unnecessarily restrict the time step because the fol- 
lowing time step is always of a magnitude commensurate with the new defined rates of 
stress and strain. 

It is of interest to note that this change primarily originated from problems pre- 
viously encountered in gas equilibrium calculations which were subsequently alleviated 
by the modifications made to subroutine SGAS. However, this change was retained 
because of potential savings in computer time. 

Self-initiated time step reduction . - The program often requires a reduction in the 
time step as previously discussed, and can require several iterations to find a-useable 
one. Consequently, a technique was employed by which the convergence to a useable 
time step is greatly accelerated. This technique is a greatly modified method of suc- 
cessive approximations, and occurs in subroutine STIM. The procedure is as follows: 
The time increment or time step called DTIME is the actual time step being used . All 
the appropriate calculations are made using this time step, including the rates of stress 
and strain. The inputed allowable stress and strain changes are then divided by these 
stress and strain rates, respectively, to obtain a new time step called DELTA. This 
DELTA represents the time step which could be used in order for the stress and strain 
changes to just equal the allowable . If the DTIME used is greater than DELTA, which 
implies that the stress and strain changes would have exceeded the allowable, DTIME is 
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set equal to DELTA and all the calculations are redone. Both the time step DTIME and 
the required time step DELTA are stored for any two consecutive iterations. If DTIME 
is less than DELTA, the stress and strain changes will be less than the allowable. In 
this case DTIME is kept as a reasonable time step and DELTA is used for the next time 
step. It is important to note that the DELTA does not remain constant, which is the 
reason that more than two iterations are required. The DELTA is dependent upon 
DTIME in a nonlinear fashion. However, the optimum time step is one in which the 
stress or strain change just equals the allowable change. Therefore, the optimum 
DTIME is one which is just equal to DELTA, and this is the basis for which a new esti- 
mated time step is determined. In other words, when two consecutive DTIME 's and the 
resulting DELTA'S (since DELTA is actually indirectly dependent on DTIME) are known, 
a slope can be determined and extrapolated to a condition where these parameters are 

equal. This technique has greatly reduced the number of iterations and reduced com- 

•¥ 

puter time . 

Allowable stress and strain rates . - The stress and strain rate during a given time 
step must be less than, or equal to, that value which, if multiplied by the time step, 
would result in an incremental change equal to the allowable values specified in the pro- 
gram input data. For example, it may be decided that it would be unreasonable to allow 

p 

a step change in the stress of more than 1000 psi (689 N/ cm ) or a strain change of more 

-3 -3 

than 10 inch per inch (10 cm/cm). 

Furthermore, the stress and strain changes in the CYGRO-2 computer program 
during a given time step are each composed of two parts. These two parts, as defined 
in reference 3, can be expressed as parts of the total stress and strain change as 

AS = AS. + S At 

Ae = Aej + e At 

respectively, for the stress change and strain change; and where ASj and Ae^ are 
stress and strain changes assumed to take place at the beginning of the time step (ref. 3, 
p. 23), and the stress and strain changes S At and e At are those taking place during 
the time step. (All symbols are defined in appendix A, ) 

Originally, the program compared only the second terms, S At and e At, with the 
allowable stress and strain changes. However, a debug of the program revealed that, 
the second terms are often equal in magnitude but sometimes opposite in sign to the first 
terms. Because of this, the time step can be unnecessarily restricted in some cases or 
allow too large a total change AS or Ae in other cases. Consequently, the program 
was modified to compare the total stress and strain change with the allowable value. 

This is accomplished in subroutine STIM. 
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Special restrictions on change in equivalent stress . - The subroutine STIM contains 
certain sign restrictions on the change in equivalent stress. For example, if the pro- 
jected equivalent stress is negative, with certain exceptions, the time step may be re- 
duced. The technique by which the time step is reduced can be relatively unstable when 
the equivalent stress is near zero; such as at the beginning of the problem or when a 
particular node (or ring) is near zero equivalent stress. It was possible to reduce the 
time step by orders of magnitude under such conditions. This was corrected by letting 
the variable containing the initial equivalent stress be the greater of either its natural 
calculated value or 500 psi (345 N/cm ). Consequently, the time step is not reduced 
when the change in equivalent stress and the equivalent stress are small. However, 
this procedure still may result in an unnaturally small time step when by chance the 
variable is negative but slightly greater than 500 psi (345 N/cm ) in magnitude. It is 
not known just how often this may occur, and further modification may be required. 

. Time step at closure of time interval between time cards . - As the problem time 
approaches that of the value on the input time card, a check is made to see if the re- 
maining time is less than the smallest allowed (TINY). If so, the time step DTIME is 
set equal to TINY and is used for the remainder of the calculations. Consequently, the 
final time can exceed the time-card time by some small amount. Since the value of the 
smallest time step allowed has been increased, as previously discussed, this section 
was slightly modified to use, in this case, either TINY or the remaining time, which- 
ever is smaller. This modification also occurs in subroutine STIM. 

Calculation of strain rates . - The total generalized strain £g is obtained by inte- 
gration of the generalized strain rate €g. The value of is approximated from the 
initial values of the independent variables and from the rates of change of these vari- 
ables. The calculated strain rate is then assumed to be constant during the time step. 
All the calculations used in determining this are in subroutine SEGT. Originally, 
subroutine SEGT contained two separate and independent methods of obtaining Only 
one ever appeared to be utilized. A logic statement always routed the calculations to 
only one of the methods. Each of these methods and the resulting modification to the 
program are discussed. 

The first method utilized the technique outlined in reference 1. The expression for 
strain rate (ref. 1) is 




AS„ 
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Derivatives are formulated into a finite difference form. The derivative 3f 00 ./3S is 

evaluated by calculating the generalized strain rate f associated with two different 

eg 

stress levels and keeping the equivalent plastic strain constant. Because this strain 
defines the intersection of the curves for creep and for plastic flow (fig. 1(b)), this in- 
tersection also remains constant. The two stress levels are determined by taking equal 
stress increments of ±AS on either side of the equivalent stress S^, as shown in fig- 
ure 1. The magnitude of AS is the same as the maximum allowable stress change 
specified in the program input. The derivative is then formed or approximated by the 
change in f with respect to the overall change in S . The derivative 3f /3e, is 
evaluated by keeping constant and calculating the strain rate associated with a new 
equivalent strain found by adding onto the current a Ae, as shown in figure 2 . 
This effectively translates the plastic creep line to a higher stress on the creep curve, 
where the new intersection between regular creep and plastic creep is at the plastic 
stress corresponding to a plastic strain of + Ae. The Ae is the maximum allowable 
strain change specified in the program input. As shown in figure 2(b), the change in f 

6g 

is zero unless the plastic flow line is involved. The value of f eg in the equation is 
determined from current values of e , and S at the beginning of the time step. Both 
of the derivatives in the equations have been approximated and a value for can be 
calculated. 

It was found upon initial application of CYGRQ-2 that this method could result in 
erratic stresses. Not only did the stresses change erratically with respect to time, but 
also with respect to radius. This method resulted in several mathematical problems. 
First of all, the derivative 3f /3e, must be small or, when the quantity 
1 ' (9f eg/ 3e ^ C lg At is evaluated, the result is a negative number. Consequently, the 
strain rate is projected as negative, which is not physically possible. Large values 
for the derivatives are possible when the equivalent stress is high enough that increasing 
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the strain an increment Ae (at constant S g ) moves the stress- strain loci from a region 
of plastic flow (or creep) into a region of regular creep. A second problem, and one 
which occurs more often than the previous, is that the creep rate associated with 
S + AS may be many orders of magnitude greater than that of S . This becomes more 

& 5 

obvious when S is very close to yield point. Then, when AS is added, the stress is 
above the plastic stress. Combining this with the concept that the plastic flow or creep 
curve is nearly vertical (log e^ with respect to log S^), e^ must be very large at 
Sg + AS. This could lead to an unstable situation. For example, if the strain rates in 
two adjacent nodes are considerably different, unstable stress distributions could result. 
These problems could be alleviated somewhat by making the required AS and Ae very 
small. This procedure, however, would result in much larger computer times. Con- 
sequently, the subroutine SEGT was slightly modified so that only the second method 
would be utilized. This second method gave much smoother stress distributions. 

The second method of calculating the generalized strain rate e^ utilizes the amount 
of strain which would occur if the strain rate existing at the beginning of the time step 
were applied throughout the time step. Here the basic derivatives being evaluated are 
slightly different. It is assumed that the equivalence between the product (3e /3S )AtS 

• o o o 

and (3e /3S )S can be utilized so that the derivative 3e/3S can be evaluated rather 

& . o & & & 

than 0€g/ 9Sg. Initially, all that is known at the beginning of the time step is the equi- 
valent strain e^ (or e^, depending on whether it is a cyclic problem (ref. 1)), the 
plastic stress of this particular strain, the equivalent stress S^, and the strain rates 
associated with both the equivalent stress and the plastic stress. It is initially assumed 

that the strain rate f occurring at the beginning of the time step occurs throughout 

eg 

the time step. The new value of the strain at the end of the time step is determined and 

the required stress to maintain this constant strain rate is solved for. This essentially 

moves the intersection on the creep - plastic flow (log S - log e ) curve to a higher 

& & 

stress, as shown in figure 3. The new required stress is that found on the new creep - 
plastic flow curve (fig. 3 (b)), keeping e^ constant. The intersection is found from the 
plastic stress-strain characteristic diagram. These values of and e^ are defined 
as the upper limits, with lower limits of the differences yet to be defined in order to 
evaluate the derivative 3e^/ 3Sg. Now a strain rate must be assumed which will result 
in the lower limits. A higher rate can be assumed for the purpose of finding 
however, it is more probable that a lower rate, resulting in lower strain, is more ac- 
curate because of such factors as stress relaxation. CYGRO-2 uses a lower strain rate. 
The lower strain rate is assumed to be the arithmetic average of the initial strain rate 
and the final strain rate at the end of the time step. This new lower strain rate, and 
consequently, lower strain, defines a new log S - log e intersection; and consequent- 

fc? & 

ly, there is a definite stress required to obtain this strain rate. The derivative can 

then be formed. If the difference between the upper and lower strains e_ - e_ * 

&? U P g, low 
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is more than the Ae requirement specified in the input data, the strain e^ is set 
equal to the upper strain minus Ae. The appropriate strain rate is determined which 
will satisfy this new lower strain. If the stress difference - S g violates AS 

requirement specified in the input data, the new lower strain is modified by another fac- 
tor to move the lower strain nearer to the upper strain. After both of these require- 
ments are met, the derivative 3 e / 3S is formed. By using the difference between 

& & 

these upper and lower limits of strain, the equivalence previously discussed is satis- 
fied; that is, 

ffg At _ € g> UP " € g> low At _ € g, up At " £ g, low At _ € g, up " € g, low _ 
aS g S g, up " S g, low S g, up " S g, low S g, up " S g, low 9S g 

where the difference between the upper and lower strain limits is used, rather than a 
difference from the strain at the beginning of the time step. This second method uses 
essentially the equation 


'g 


= A 


gt 


A S 
gg £ 


where 



= € 


g 


S 


g 



0 € 

as 


1 

g 


and where de^/ is evaluated as previously described, and A^^. is evaluated by 
holding S constant and determining what the effective strain rate must be. This is 
not simply a constant creep rate, especially if the stress is in the plastic flow region. 
(This is because of the shifting of the intersection on the log - log k^ curve. ) Once 
the strain is determined for constant S , the strain rate may be determined. The 
strain is determined by utilizing the preceding derivative for A^^. The strain is ex- 
trapolated backwards to from the upper stress using the slope. The strain rate is 
then determined from this strain. 

It was found that this procedure always underestimates the amount of strain rate 
e and, in fact, sometimes extrapolates to negative strain rates. This appears in the 
CYGRO-2 output as zero strain rate because only strains greater than, or equal to, 
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zero are allowed. A check revealed that this can be very seriously in error because the 
curve is assumed to be linear for extrapolation. However, an evaluation shows it to be 
extremely nonlinear; consequently, when the extrapolation occurs over large stress 
changes, a large error occurs. This is shown in the following sketch, where the solid 
line represents the extrapolation and the large-dashed line represents the true behavior: 



Therefore, the subroutine was modified to obtain more accurate estimates for the strain 
rates. It was felt that the easiest way to attack this problem would be to modify the sec- 
ond method, since the only changes required would be to use an alternate method of de- 
termining the strain rate at constant S^. The first method was retained within the pro- 
gram although it is not utilized. 

The strain rate is again evaluated from the equation 


e„ = e. 


d€„ . 

+ — So- 

So- dS cr g 

g g 


where the value of the strain rate assuming constant stress can be obtained easily from 
the creep curves. An examination of the following sketch shows how this can be done by 
an iteration technique: 
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It must be remembered that the intersection between the regular creep line and the ac- 
celerated creep or plastic flow line is affected by the creep assumed to have taken place 
during the time step. Consequently, if the stress S is above the plastic stress, an 

w 

iteration scheme is required; if below, no iteration is required. The object is to find 
the value of strain rate (indicated by k) while keeping constant, which satisfies both 
the constant stress and change in position of the intersection due to the change in strain. 
It is obvious that only one such solution exists . An iteration scheme finds this solution 
in just a few iterations . A value for the first term in the strain rate equation is then 
known. 

Next, the value of the derivative de^/ dS^ must be found. This procedure was left 
nearly identical to the original CYGRO-2, with the exception that the lower value of the 
strain rate becomes that which was just solved for. Consequently, this slope is based 
upon strain rates which are closer to those just deviating from the condition of constant 
stress S„. 

With this scheme there are no resulting zero strain rates with finite stresses, and 
the slight increase in iterations in this small area adds no significant amount of compu- 
ter time. It might be mentioned that previous to this change a particular situation 
existed such that the stress in one of the nodes of a sample problem was increasing but 
the creep rate was zero. After a few time steps the stress became so large an overflow 
occurred. 


Optional Modifications 

The CYGRO-2 program was further modified to give additional calculations and 
additional printout of the calculated parameters than what was done in its original form. 
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It was also modified slightly to give a different comment card format for the data deck 
and a slightly different printout format. 

Input titles for the data deck were altered. Any number of comment cards can be 
read in and printed previous to the data. The number of title cards is controlled by the 
first two columns on each card. As long as the space is left blank, continual reading of 
additional comments cards occurs. The last card then contains a number and the next 
card must be a regular data card (card 900). This allows a very descriptive set of com- 
ments so that different problems being run are easily distinguishable. 

The printed output which was divided into groups was modified in two ways. The 
write format was changed so that the data are spread out more; in other words, simply 
more space between numbers. Secondly, more frequent printout is obtained. Originally, 
a printout of groups 3 to 7 (ref. 1) occurred at each time card. Now a counter called 
MMM controls additional printout between time cards of group 4, and a counter called 
NNN controls the printout between time cards of groups 4 to 7. These are essentially 
print frequency counters. They are specified on one of the input cards: card 907, items 
7 and 8. For example, if MMM equals 5 and NNN equals 10, every fifth time step, 
regardless of elapsed time, would result in group 4 being printed and every 10th time 
step, groups 4 to 7. The program still prints the entire group (3 to 7) when the elapsed 
time equals that specified on the time cards. The additional printout allows the calcula- 
tions and progress to be followed much better. Group 4 contains time, critical dimen- 
sions, maximum stresses and strains, and contact forces. Groups 5 to 7 include indivi- 
dual ring data on temperatures, radii, stress and strains, and volume growths, and 
bubble information. 

One additional printout of all data occurs at the point in time where the fuel and clad 
have just contacted. These data provide useful information on free fuel growth and exact 
contact time. 

An additional output parameter was included in the group 4 printout which indicates 
whether plastic flow (accelerated creep) or regular creep has taken place in each of the 
rings of fuel or clad during the time step. The plastic flow index is zero for creep, a 
normal condition until the equivalent stress exceeds the plastic yield stress S . Then 

r 

the index is equal to 1. 

Finally, a fractional fuel volume change and fractional radial dimensional change of 
fuel and clad were added to the group 4 output. These are based on instantaneous dimen- 
sions relative to the original dimensions at time zero. The diametral changes are cal- 
culated for the fuel inner diameter and for the fuel and clad outer diameters. Volume 
changes are calculated for the fuel only. 
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Application to Lewis IBM 7094 


Storage problems . - Initially, the CYGRQ-2 program exceeded the 32 000-word core 
storage of the IBM 7094. This storage problem was somewhat alleviated by eliminating 
reference to any of the Westinghouse- Bettis environmental subroutines which were avail- 
able on the Lewis computer. Elimination of these subroutines would have been done 
regardless of a core storage problem. Two of the Bettis environmental subroutines had 
no corresponding Lewis routine. Regular subroutines were made of these (TITLES and 
LINES) and added to the CYGRO-2 deck. In addition, the subroutine calling the compu- 
ter time (TICK) was modified in name to call the Lewis timer (TIME1). All the call 
statements to eliminated subroutines were removed. The subroutine LATON was also 
split into two subroutines to help decrease required storage. These were called LATON 
and LATON1 . The subroutine LATON1 is simply a section of LATON which allowed an 
alternate method of reading in the operating history. This resulted in the CYGRO-2 
program having a main program and 20 subroutines. To further reduce the required 
core storage to that amenable to the Lewis computer, four of the subscript ranges of 
variables used in data storage for restart problems were reduced (FORTRAN names 
MUF1, MUF2, BUF1, BUF2). These are no longer used for anything in the modified 
program and could probably be removed if necessary. 

The CYGRO-2 program fit onto the computer with only a small amount of storage to 
spare. It was later found that when further modifications to the program were necessary, 
core storage continually became a problem. Hence, it was decided to eliminate addi- 
tional subroutines dealing with tape storage of data because the advantage gained in core 
storage would be more beneficial, especially since a modified printout option allowed 
more frequent printout of the parameters. There are currently 13 subroutines in addi- 
tion to the main program in CYGRO-2. The subroutines are SEND, SEGT, SRAW, STIM, 
CTAC, SGAS, STEM, SEAS, FINIS, SKEN, LUST, LATON, and LATON1. 

General FORTRAN modifications . - A few general modifications were required to 
make the FORTRAN compatible with the Lewis computer. One of these was the call to 
the computer timer previously mentioned. Others included redundant parentheses, in- 
correct reference to read-and- write tapes, incompatible data statements, and an incom- 
patible field width for the A- conversion format on title reading and writing. 

Overflows and underflows . - Several modifications were required to eliminate over- 
flow and underflow conditions on the IBM 7094. These always occurred in, or arose 
from, calculations performed in the subroutine SEGT. 

First, the creep relations and the stress-strain relations in subroutine SEGT were 
originally written in the form 


e = exp 


A lnlO + B In 
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where e is the creep rate or the plastic strain for a given stress S. This equation was 
modified in form to that originally given in reference 1, or 


e = 



This eliminated underflows that normally occurred when the exponential was evaluated. 
The second form circumvents this problem by evaluating two lesser exponentials and 
then directly multiplying them. 

Secondly, many of the parameters involved in strain and strain rate calculations 
were double precisioned. Double precision avoids what would naturally be a computer 
zero when very small differences are taken between large numbers. These parameters 
were mostly strain parameters in the subroutine SEGT. 

Finally, there are a few critical equations in SEGT where two or more small num- 
bers are calculated and then multiplied together, resulting in an underflow. This is 

particularly true for very small time steps. This problem can be corrected by checking 

-38 

the magnitude of the variables involved; and if an underflow would occur (<10 ), the 

variable is set equal to zero. (An appropriate control card could accomplish the same 

thing but only a set number of times before the problem is terminated. ) Low values of 

-38 

strain or strain rates, less than 10 , can be set equal to zero with no ill effect. 


SAMPLE PROBLEM 
Description of Problem 

The CYGRO-2 computer program is being used to predict fuel swelling for uranium 
mononitride fuel pins. Therefore, the sample problem included here is a UN-fueled fuel 
pin. The geometry of the fuel pin is as shown in figure 4. It consists of a cylindrical 
fuel pin with UN fuel and T-lll (Ta-8W-2Hf) clad. The pin is hollow, with a fuel-clad 
gap of 1. 5 mils (3. 8x10 cm) on the radius. The UN fuel is hot pressed and sintered to 
a density of 96 percent. Power and fuel burnup requirements are 9172 Btu per hour per 

O 

cubic inch (164 W/cm ) and 2.77 percent burnup on 50 000 hours. Burnup here refers to 
fissions only and not total fuel atoms destroyed. An additional heating power of 1 percent 

O O 

is assumed within the clad material (91. 72 Btu/ (hr) (in. ); 1. 64 W/cm ). Temperature 
of the clad surface is assumed to be 1720° F (1211 K). 

The gas-gap gas is assumed to be helium. The thermal conductivity of helium at 
1720° F (1211 K) is approximately 0. 0183 Btu per hour per inch per °F (0. 0038 
W/ (cm) (K)) (ref. 9). It was felt that although the xenon fission product would contribute 
heavily to the gas within the void, the fuel- clad gap would close before an appreciable 
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error in fuel temperature would occur. It has also been assumed that the interface 
thermal resistance upon fuel-clad contact is negligible. Therefore, these contact con- 
ductivities are very high. 

Experimental evidence supports the assumption that only about 5 to 10 percent of the 
fission products which are gaseous escape from the UN fuel (ref. 10). For this problem 
5 percent was assumed. Only one type of bubble has been assumed, although up to five 
types are possible in the CYGRO-2 program, and the bubble density was assumed to be 
of the order of 2.7x10*^ bubbles per cubic inch (1. 6xl0 15 bubbles per cm^). Experi- 
mental data of various fuel specimens of a variety of composition has revealed an ob- 

13 

servable bubble density as low as 5. 5x10 bubbles per cubic centimeter to as high as 
17 

10 bubbles per cubic centimeter (refs. 11 to 13). Theoretical nucleation densities for 

IB 1 ft 

inert gas vary between 10 to 10 sites per cubic centimeter (refs. 14 and 15). 

Van der Waal’s equation-of-state exponents for xenon were used for the gas bubbles. 
Surface tension on the gas bubble - fuel interface was assumed to be 500 dynes per cen- 
timeter (2. 855x10 lb/in. ). Gas generation within the bubble is from two sources. 

First, there is a gaseous fission product generation which is temperature dependent. 

For this sample problem, only the stable and inert gases were considered. Approxi- 
mately 15 percent of the fission products are gaseous (0. 3 molecule of gas generated for 
each fission event). Sources for this data are references 16 and 17. A detailed study of 
resulting fission products is given in reference 18. Most of the nitrogen either recom- 
bines with the fuel or combines with other fission products. A small fraction is left un- 
combined, however, and adds to the gas bubbles. This amount (ref. 19) is temperature 
dependent but for this sample problem is assumed to be 0. 09 molecule per fission event. 
The total gas generation rate is, therefore, 0.39 molecule per fission. 

Solid fission product expansion was assumed to be 0. 85 percent solid volume growth 
per percent burnup. This number has been nominally calculated on a theoretical basis 
for UQg as 0. 864 percent growth per percent burnup (ref. 20). Based on 0.3 molecule 
of fission gas produced per fission and another 1. 0 molecule of fission products per fis- 
sion being used to fill the voids of uranium destroyed, the increase can be shown to be 
0. 7 percent swelling per percent burnup based only on physical space requirements. Of 
course, the fission products will not exactly fill the voids. The value chosen for this 
example problem was 0. 85 percent AV/V per percent burnup. Isotropic fuel expansion 
was assumed throughout. 

All the mechanical properties were considered to be independent of fission flux or 
flux rate. Both the mechanical and appropriate physical properties were assumed to be 
temperature dependent, with the exception of the Poisson’s ratio. Poisson’s ratio for 
uranium mononitride is 0.28. (A value of 0. 26 from ref. 21 was revised to a more dense 
material. ) No data exist for Poisson’s ratio for T-1X1 alloy. Therefore, a ratio of 0.35 
was assumed because the alloy should not deviate much from the parent material tantalum 
(ref. 22). Other properties used for T-lll are modulus of elasticity (refs. 23 and 24), 
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coefficient of thermal expansion (refs. 23 to 26), and the thermal conductivity (refs. 23 
and 27). These properties are given in table I. Uranium mononitride properties as used 
in the program are given in table II. Only one data point for the modulus of UN was 
available (ref. 28). Temperature dependency was estimated based on the limited amount 
of known temperature dependency of uranium dioxide and uranium carbide as reported in 
reference 28. Thermal expansion data are available in references 28 and 29. The 
source of data on thermal conductivity of UN is reference 30. 

The form of equation which describes the stress-against-strain characteristics of a 
material is 


e 


P 



where 

e plastic strain, in. /in. (cm/cm) 

r 

o 

S stress at plastic strain e , psi (N/cm ) 

r r 

P, Q characteristic exponents, dependent on temperature and flux 
Similarly, the form of equation describing creep behavior is 


‘g 



where 

e equivalent strain rate, in. /(in. ) (hr) (cm/ (cm) (hr)) 

o 

o 

Sg equivalent stress, psi (N/cm ) 

C, G characteristic exponents, dependent on temperature and flux 

The stress-against-strain characteristics of T-lll are shown in figure 5. These curves 
are the result of both published data (ref. 31) and unpublished TRW data (private com- 
munication from K. D. Sheffler, TRW, Inc., Cleveland, Ohio). Since the form of the 
preceding equation only approximates the true curve, only the characteristic curves and 
the 0.2 percent yield points are shown in figure 5. The corresponding values of P and 
Q are given in table II. Creep properties for T-lll were obtained from reference 32. 
The exponents are tabulated in table III. 

There appears to be no data available for the strength of uranium mononitride. It 
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was assumed that since there was some high-temperature creep data available at 

O 

8000 psi (5520 N/cm ), the yield point must be at least above this. The assumed char- 
acteristics and exponents are given in figure 6 and table IV, respectively. Creep 
characteristics are similar to those in reference 33 . The creep equation exponents for 
UN are given in table IV. 

It is not intended that the property values used in the sample problem be interpreted 
as correct or up to date; they are only the specific values used. 

Only one power cycle is considered over a period of 50 000 hours. The tempera- 
tures and power come up to maximum in 1 hour and remain there for the duration of the 
problem, then ramp down to zero power and room temperature in 1 hour. It is assumed 
that the coolant pressure doubles upon heating and that the gap pressure increases with 
absolute temperature and ramps linearly while hot from 73. 5 to 368. 0 psi (50. 7 to 
253.7 N/cm 2 ). 


Input Data Card Description 

A considerable amount of experience has been gained in both debug computations for 
making modifications to CYGRO-2 and in making studies correlating experimental data 
with CYGRO-2 predictions (ref. 34). It may be of benefit, therefore, to discuss the input 
data and the resulting potential problems encountered in using CYGRO-2. 

In general, the input is necessarily lengthy and complicated. Because of its length 
it may take a lengthy period of time to obtain and put the data into useable form. The 
program user must not only obtain the material's mechanical and physical properties, 
specify the fuel-pin geometry, and determine the time -dependent conditions, but must 
correctly model the problem (such as the number and distribution of bubbles) through the 
many modeling options available . This can be a disadvantage in that the user must be 
knowledgeable about fuel swelling in order to correctly choose conditions. And, if a new 
material is being considered, such as some refractory alloy, the required properties are 
not always known. This further complicates the problem. 

Further explanation and comments of input data cards, as they appear in the modi- 
fied CYGRO-2 program, are given in table V. The program input is given in appendix B. 


Sample Problem Results 

A FORTRAN statement form of the input data deck is given in appendix B. This par- 
ticular sample problem has a total of 89 data cards, including a blank end- of- problem 
card. Also included in appendix B is a sample problem computer printout. This print- 
out contains (1) the identification page; (2) the input data, including both the comment 
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cards and data cards; (3) a program -calculated table of the plastic stress against strain 
and creep stress against strain rate; (4) all the calculated parameters at full power 
(1 hr); (5) an abbreviated printout at the time of fuel-clad contact (18 907 hr); and (6) a 
printout of all the parameters at the end of life (50 001 hr). 

Total fuel swelling for this problem was 5.2 percent, of which 2.35 percent is solid 
fission product swelling. There is a range in local swelling, ranging from 6. 0 percent at 
the inside node of the fuel to 4. 6 percent at the outer node. There is a considerable dif- 
ference in movement of the fuel between nodes, however. The inside surface of the fuel 

-2 

has been pushed into the central void by about 10 mils (2, 54X10 cm), whereas the outer 

-3 

surface moved out only 1 mil (2. 54x10 cm). For this particular pin the central void 
serves as a sink for swelling, resulting in very little clad deformation - an ideal 
situation. 

The sample problem ran in a computer time of approximately 10 minutes. Problems 
involving variations in fuel loading, burnups, irradiation times, cycling, fuel-pin diam- 
eters and gaps, and solid or hollow pins have taken from as little as 2 minutes computer 
time to as high as 60 minutes. Problems which seem to require additional computer 
times are those involving high gas generation rates and high temperatures. This sample 
problem used 833 time steps for 50 001 hours. 

It is characteristic of the problems run thus far to only load the clad enough to cause 
plastic flow during a thermal expansion. During the rest of the time the clad is in creep. 
Running the same problem using different sets of time cards can result in variations in 
stress when the total strains are the same. This is believed to be due to the stress levels 
being very sensitive to very small changes in strain and is not significant in determining 
fuel-pin swelling. 


CONCLUDING REMARKS 

The CYGRO-2 computer program as initially received was not compatible with the 
Lewis IBM 7094 computer. The program was subsequently modified such that it would be 
compatible. It was intended to use the program basically as a design tool after compar- 
ing theoretical and experimental fuel swelling. Many computer problems were encoun- 
tered with the program. Consequently, a debugging and program modification procedure 
was instituted. Some of these modifications alleviated specific calculational problems, 
whereas others improved accuracy and computation time. 

A discussion of the input data serves to clarify some of the variables that were of 
questionable nature and magnitude when the program was first utilized. It was found that 
the tabular functions can, if not carefully set up, result in error and, in some instances, 
prevent completion of the problem. 
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A sample problem of a UN fuel pin, with a detailed discussion of the input data, has 
been included. This sample problem has been utilized to show the changes that have been 
made in input and output data. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, January 13, 1971, 

120-27. 
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APPENDIX A 


SYMBOLS 

coefficient exponent in strain and strain rate equations, replaces C and P 
rate of change of strain rate with respect to stress rate 
strain rate at constant stress 

stress exponent in strain and strain rate equations, replaces G and Q 
coefficient exponent in creep equation 

matrix relating equivalent cyclic strain to generalized equivalent strain 

strain or strain rate 

generalized strain rate function 

stress exponent in creep equation 

coefficient exponent in plastic stress- strain equation 

stress exponent in plastic stress- strain equation 

stress 

stress rate 

stress change 

generalized equivalent stress 
generalized stress rate 
plastic or flow stress 
time 

total volume 
gas bubble volume 
strain 
strain rate 
strain change 
generalized strain rate 
equivalent plastic strain 
plastic strain 



Subscripts: 

apparent 

i 

low 

true 

up 


apparent value 

change in variable assumed to take place suddenly at beginning of time step 
lower value 
true value 
upper value 
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Sample Problem Output 


STOP f irt£=30« PAGES -2 00 9 FORMS-PLAIN 

SOVRFLW CONTINUE 

SMAJUND CONTINUE 

SZERDIV CONTINUE 

$ I B JOB GO® NOSOURCE 


0 SAMPLE PROBLEM - LEWIS MODIFIED CYGRO-2 
0 TEMP CLAD SURFACE = 1720 F 
0 T-lll CLAD 

0 UN FUEL 96 PERCENT DENSE 

0 FUEL ID = 0.300 OD = 0,621 

0 3 MIL DIAMETRAL GAP 

0 PIN OD = 0,750 

0 INITIAL HOT GAS PRESSURE = 73.5 PSI - FINAL = 368. PSI 
0 2,77 PERCENT BURNUP 50000 HOURS 

0 BUBBLE DENSITY 2.7E+16 

0 OXYGEN CONTENT OF FUEL 400-500 PPM 

0 5 PERCENT GAS LEAKAGE 

0 0,85 PERCENT SOLID GROWTH PER PERCENT BURNUP 


3 


wcs 

FID 

FI 2) 

FI 3 ) 

FOOT ( 1 ) 

FD0TI2) 

FOOT 1 3 ) 

0 

9.00 + 2 

3.0000E+01 

0. 

1.0000E+01 

UOOO0E+O3 

UQ0OOE-O6 

1.0000E-04 

1.0000E-02 

5.3000E-01 

9.0000E+02 

TG(K ) 

DFTIDK) 

DFT I 2 ®K ) 

DFTI 3*K) 

RGASNIK) 

ZUYKB 

0 

0 

9« OtC + 2 

0. 

1 « OOOOE+OO 

1 » OOOOE+OO 

1, OOOOE+OO 

1. OOOOE+OO 

0, 

0. 

0. 

9.010DE+02 

1.0000E+03 

1 .OOOOE+OO 

1. OOOOE+OO 

1. OOOOE+OO 

1. OOOOE+OO 

0. 

0. 

Oo 

9.0200E+02 

3, OOOOE+03 

1 e OOOOE+OO 

1. OOOOE+OO 

1 » OOOOE + OO 

l. OOOOE+OO 

0. 

0. 

0. 

9.0300E +02 

CGAP 

HGAP 

H NOGAP 

HWATER 

RIVDW) 

AIVDW) 

8IVDW) 

TABS 

9.04 + 2 

1.8300E-02 

1 « 0000E+10 

1 » 0000E+10 

1.0000E+10 

1 « 8540E+04 

1.4200E+Q3 

4» 6300E+10 

4.6000E+02 

9. 0400E+02 

VAN DER WAALS LOWER TEMPERATURE LIMIT IS -20. 




- 


FCORE 

FCONIR) 

FCONIZ) 

DS I MAX ) 

DEIMAX) 

RPLUG 

QPLUG 

TPLUG 

9.05 + 2 

-0* 

-0. 

-0. 

2. OOOOE+03 

1.0000E-03 

-0. 

-0. 

5 © OOOOE+O 3 

9.0500E+02 

SGMIN 

MPQ 

MCORE 

MGAS 

MGL 

EOOGS 

MFILE 

PFILE 

9® 06 + 2 

1. OOOOE+OO 

0. 

. 1 « OOOOE+OO 

0 . 

1. OOOOE+OO 

0 , 

0 . 

1. OOOOE+OO 

9® 0600E+02 


-1 


MCEE 

DVMAX 

HIST 

RAP 

RBP 

MAXC 

mn 

NNN 

9® 07 + 2 

0* 

1 0 OOOOE-03 

0® 

l.OOOOE+OO 

1® 3000E+00 

0® 

UOOOOE + 03 

1 * OOOOE+O 3 

9® 0700E+0 2 


MEAT 

PRQPERT IES- 

— 






NR 

NG 

RSC 

RSP 

MCGPQ 

OGMAX 

EZ 

S»TEN® 

I .000 + 3 

1® 1000E+01 

1 * OOOOE+OO 

1 • OOOOE+OO 

1 ® OOOOE+OO 

2*. OOOOE + OO 

9.1720E+D3 

0., 

2® 8550E-0 3 

1, OOOOE+03 
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RADIUS ( KR) 

SR(KE) 

SC(KE) 

SZtKE) 

EGC(KE) 

EGZ(KE) 

EG(KE) 

EVtKE) 

I. 001 + 3 

1* 5000E-01 

-0® 

-0, 

-Oa 

-0. 

-0® 

-Oa 

-0® 

l s 0010E+03 

lo 5300E-01 

-0* 

-0* 

-0® 

-Oa 

-Oa 

-Oa 

-0® 

1 ® 0010E +0 3 

1« 5500E-01 

-0. 

-0® 

— 0 e 

-Oa 

-Oa 

-Oa 

-0® 

la 0010E+03 

1. 6050E— 01 

-0® 

-0, 

-Oa 

-0® 

-Oa 

-0 * 

-0* 

1 * 0010E + 03 

U8190E-01 

-0* 

-Oa 

-0* 

-Oa 

— 0 a 

-Oa 

-Oa 

1.0010E+03 

2» 0330E-01 

-0® 

“0 « 

-Oa 

-0® 

-Oa 

-Oa 

-0. 

1*00 10E+0 3 

2 * 2480E-0 1 

-0® 

-0 « 

“Oa 

-Oa 

-0* 

-Oa 

-Oa 

1» 0010E + 03 

2*46206-01 

-0® 

-0. 

-Oa 

— 0 a 

-Oa 

-0® 

-0* 

1 aOOlOE + O 3 

2* 6760E-01 

-0® 

-Oa 

-Oa 

-0. 

-Oa 

-Oa 

-0* 

1.0010E+03 

2® 8900E-01 

-0* 

-Oa 

— 0 o 

-Oa 

-Oa 

-Oa 

-Oa 

1.0010E+03 

3* 1050E-01 

-0* 

-Oa 

-Oa 

-0® 

-Oa 

-o « 

-0® 

1 ® 0010E + 03 

MOLES ( KG ) 

AM I KG ) 

BM(KG) 

VH(KG) 

V ( KG) 

NH(KG) 

NNHTtKG) 

NHT ( KG) 

I « 002 + 3 

6® 4500E-08 

1 • 426QE-07 

2 a 0000E-06 

4* 0000E-D2 

1 » OOOOE+OO 

2 o 7000E + 1 6 

— Oa 

2. 7000E+1 6 

1* 0020E + 03 

RALPHAtR) 

RALPHA(C) 

RALPHAtZ) 

RV(R) 

RV{ C) 

RVtZ) 

0 

0 

I e 003 + 3 

1 * OOOOE + OO 

1 ® OOOOE+OO 

l.OOOOE+OO 

3* 3300E-0 1 

3 • 33006-01 

3 e 3300E-0 1 

-0® 

-0® 

U0030E + 03 

Ptl) 

P{ 2) 

P( 3) 

0 

0 

0 

0 

0 

1.004 + .3 

-0® 

6® O000E+03 

1 a GOOOE+04 

-Oa 

-Oa 

-0. 

-0. 

-Oa 

1® 0040E+03 

kr(=u 

RQ ( 1) 

KR 

RQ(KR) 

KR 

RQ t KR ) 

KR 

RQtKR) 

1.005 + 3 

1® OOOOE+OO 

1 « OOOOE+OO 

3 a OOOOE+OO 

1® OOOOE+OO 

6« OOOOE+OO 

1 » OOOOE+OO 

lalOOOE+Ol 

1. OOOOE+OO 

1.0050E+03 

TEMPT ( J ) 

AL PHAT ( J ) 

ET(J) 

NUT t J ) 

CTU) 

RVP ( 1 » J ) 

RVP ( 2 » J ) 

RVP { 3* J ) 

I o 0 J 0 + 3 

0 

CFT t K « J ) 

GFT { K » J ) 

PFT ( K» J ) 

QFT ( K* J ) 

0 

0 

0 

IaOJK + 3 

3® 2000E+01 

4« 670QE-06 

3 « 1000E+07 

2« 8000E-01 

7® 7000E-01 

1« OOOOE + OO 

1 ® OOOOE+OO 

1 « OOOOE+OO 

1 a 0100E+03 

-0® 

-3a OOOOE+O 1 

5 » 9020E+00 

-2« 2699E+01 

2« OOOOE+Ol 

-0® 

-0. 

-0. 

1 .0 110E+03 

-0® 

-3* OOOOE+O l 

5 ® 9020E +00 

-2. 2699E+01 

2 «OOOOE+Ol 

-0. 

-0. 

-Oa 

1 ®0 120E+03 

-0® 

-3 « OOOOE+Ol 

5 • 9020E+00 

-2. 2699E+01 

2. OOOOE+Ol 

-Oa 

-0. 

-0. 

1 .01306+0 3 

9. 3200E+02 

4.6700E-06 

3 ® 030QE+07 

2«8000E-01 

9.8200E-01 

1. OOOOE+OO 

1 a OOOOE+OO 

1. OOOOE+OO 

1® 0200E + 03 

-0. 

-1« 5131E+01 

5 a 9020E+00 

-2.2053E+01 

2 . OOOOE+O 1 

-0. 

-Oa 

-0. 

1.0210E+03 

-0® 

-1* 5131E+01 

5 « 9020E+00 

-2.2 053E+01 

2* OOOOE+Ol 

-Oa 

-Oa 

-0. 

1® 0220E + 0 3 

-0® 

-1 « 5131E+01 

5 « 9020E+00 

— 2« 2053E + 01 

2. OOOOE+Ol 

-0® 

-Oa 

-0. 

1.0230E+03 

1® 8320E+03 

5® 0600E-06 

2 ® 8500E + 07 

2* 8000E-01 

1 • 1600E+00 

laOOOOE+OO 

1 » OOOOE + OO 

l.DOOOE+OO 

1. 0300E+03 

-0. 

-8« 5344E+0Q 

5 « 902QE+00 

-2. 1405E+01 

2® OOOOE+Ol 

-Oa 

-0® 

-0. 

1.0310E+03 

-0* 

-8 ® 5344E+00 

5 » 9020E+00 

—2. 1405E+0 1 

2.0000E+01 

-0. 

-0® 

-Oa 

1.0320E+03 

-0® 

-8. 5344E+00 

5 « 9020E+00 

-2.1405E+01 

2. OOOOE+Ol 

-0® 

-0. 

-Oa 

1® 0330E+03 

2* 7320E+03 

5 ® 3400E-Q6 

2 ® 2800E+07 

2# 8000E-01 

1 a 2480E+00 

1®OOOOE+O0 

1 » OOOOE+OO 

1. OOOOE+OO 

1.0400E+03 

-0® 

-5® 6583E+G0 

5 ® 9020E+00 

-2.0759E+01 

2. OOOOE+Ol 

-0® 

-0. 

-Oa 

1 a 0A10E+O3 

-0* 

-5® 6583E+00 

5.9020E+00 

-2. 0759E+0 1 

2 oOOOOE+Ol 

-Oa 

-0. 

-Oa 

1. 0420E+03 

-0* 

-5® 6583E+00 

5 ® 9020E+00 

-2« 0759E+01 

2 . OOOOE+Ol 

-0® 

-0® 

-0® 

1® 0430E + 0 3 


NR 

CLAD 

NG 

PROPERTIES- 

RSC 

RSP 

MCGPQ 

QGMAX 

EZ 

SaTENa 

I a 0 00 + 3 

6.0000E+00 

0. 

1. OOOOE+OO 

1 . OOOOE+OO 

2 « OOOOE+OO 

9 a 1 720E + 0 1 

Oa 

2® 8 550E-0 3 

2. OOOOE+O 3 

RADIUS(KR) 

SRtKE) 

SCtKE) 

SZtKE) 

EGCtKE) 

EGZIKE) 

EGtKE) 

EVtKE) 

I. 001 + 3 

3® 1200E-01 

-0® 

-0. 

— 0 a 

-0. 

-0. 

-Oa 

-0® 

2 a0010E + 03 

3. 1700E-01 

-0® 

-0® 

-0. 

-0. 

-Oa 

-0® 

-Oa 

2. 0010E+03 

3. 3200E-01 

-0® 

-0* 

-0. 

-0. 

-0® 

-Oa 

-Oa 

2 ® 0010E + 03 

3*47006-01 

-0® 

-0* 

-0. 

-Oa 

-0. 

-0. 

-Oa 

2.00106+03 

3® 6200E-01 

-0. 

-Oa 

-0® 

-Oa 

-0. 

— Oa 

-Oa 

2.0010E+03 

3® 7500E-01 

-Oa 

-0* 

-0® 

-Oa 

-Oa 

-Oa 

-0. 

2« 0010E + 03 

RALPHAtR) 

RALPHAtC) 

RALPHA(Z) 

RVtR) 

RV ( C } 

RV(Z) 

0 

0 

1.003 + 3 

1 ® OOOOE+OO 

1. OOOOE+OO 

1 . OOOOE+OO 

3® 3300E-01 

3® 3300E-01 

3* 3300E-01 

0* 

Oa 

2® 0030E + 0 3 

Ptl) 

P(2) 

P t 3 ) 

0 

0 

0 

0 

0 

1.004 + 3 

0® 

6* OOOOE+03 

1 « 0000E + 04 

— 0 a 

-0. 

-Oa 

-0 a 

-Oa 

2« 0040E + 03 

KR { = 1 ) 

RQ ( 1) 

KR 

RQ(KR) 

KR 

RQ(KR) 

KR 

RQ t KR) 

1.005 + 3 

1 ® OOOOE + OO 

1 * OOOOE+OO 

2 ® OOOOE+OO 

1 * OOOOE+OO 

3® OOOOE+OO 

1 » OOOOE+OO 

6. OOOOE+OO 

1 * OOOOE+OO 

2.0050E+03 

TEMPT (J) 

ALPHATt J) 

ETU) 

NUTt J) 

CTt J) 

RVP { 1 ? J ) 

RVP ( 2 ? J ) 

RVP ( 3« J ) 

I.OJO + 3 

0 

CFTtK,J) 

GFT(K»J) 

PFT ( K » J ) 

QFT t K® J ) 

0 

0 

0 

IaOJK + 3 

3* 2000E+01 

3.0000E-06 

2 * 6000E+07 

3® 5000E-01 

2« OOOOE+OO 

1* OOOOE+OO 

1* OOOOE+OO 

1 » OOOOE+OO 

2« 0 100E+0 3 

-0® 

-3® OOOOE+Ol 

3 * 6920E + 00 

-3* 1 270E+01 

1 . 5 10QE+01 

-0. 

— Oa 

-0. 

2 » 01 10E+0 3 

-Oa 

-3® OOOOE+Ol 

3 * 6920E+00 

-3c 1270E+0 1 

1.5 100E+01 

-0 * 

-Oa 

-0. 

2* 0120E + 03 

-0* 

-3 * OOOOE+O 1 

3 « 692 OE+OO 

-3* 1270E+01 

1 « 5 10QE+01 

-0, 

“Oa 

-Oa 

2 * 0 130E + 03 

1® 8320E+03 

3* 6500E-06 

2 « 3000E+07 

3»5000E-01 

2. 8100E+00 

laOOOOE+OO 

1* OOOOE+OO 

1 .OOOOE+O 0 

2.0200E+03 

-0® 

-K0474E+01 

3 « 6920E+00 

-2 * 5427E+01 

1.5100E+01 

-0® 

-Oa 

-0* 

2 .0210E+0 3 

-0* 

-U0474E+01 

3 ® 6920E + 00 

-2® 5427E+01 

1 « 5 100E+01 

-Oa 

-0® 

-0® 

2.0220E+03 

-0. 

-1* 0474E+01 

3 o 6920E + 00 

-2a 5427E + 01 

1 a 5 100E+01 

-Oa 

-Oa 

-Oa 

2*02306+0 3 
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2® 4000E+03 
-0. 

-Oa 

-0. 

3.9200E-06 
-7® 23Q0E+00 
-7 .23006+00 
-7 « 2300E+00 

2 » 2000E+07 
3» 6920E+00 
3 • 6920E+00 

3 ® 6920E+00 

3.5000E-01 
—2 a 3900E-t-01 
-2®3900E*0i 
-2.3900E+01 

3«OOOOE+O0 

1*5100E+01 

la5l00E+01 

1®5100E*01 

1 ® OOOOE+OO 

-Oa 

-0® 

-0. 

l*0000E+00 

“Oa 

“0* 

-0® 

1 • OOOOE^OO 
-0* 

-Oa 

-0* 

2.0300E+03 
2 ® 0 3 1 0 E +0 3 
2.0320E+03 
2.0330E+03 

2® 7320E+03 
-0® 

-0® 

“0® 

4® 0800E-06 
“5. 8776E+00 
-5® 8776E+00 
-5® 8776E+00 

2 • 1400E+07 
3 * 6920E+00 
3 • 6920E+00 
3« 6920E+00 

3® 5000E-01 
-2.2664E+01 
-2®2664E+01 
-2®2664E+01 

3a0500E+00 
1 « 5100E+01 
1 * 5100E-»-01 
1® 5100E+01 

1 eOOOOE+OO 

-Oa 

-0* 

1. OOOOE+DO 
-0® 

-0. 

-0® 

1 • 0000E<-00 
-0® 

-0® 

-0* 

2*0400E+03 

2.0410E+03 

2.0420E+03 

2*0430E+03 


KM KT 

1 1 

1 2 

1 3 

1 4 

2 1 

2 2 

2 3 

2 4 

KM KT 

1 1 


l 2 


1 3 


1 4 


2 1 


2 2 


2 3 


2 4 


<F 

SP(OaOOl) 

SP I Oo 002 ) 

SP ( 0.005 ) 

SP(O.OIO) 

SP10.020) 

SP < 0. 050 ) SP(O.IOO) 

SP< 0*200) 

SP (0.500 ) 

SP( 1.000) 

1 

9659® 

10000. 

10469® 

10838. 

11220. 

11746, 

12160® 

12589. 

13179. 

13644* 

2 

9659® 

10000® 

10469® 

10838® 

11220. 

11746. 

12160. 

12589® 

13179® 

13644. 

3 

9659* 

10000® 

10469® 

10838* 

11220. 

11745. 

12160® 

12589. 

13179* 

13644. 

1 

8967* 

9283. 

9718. 

10061. 

10416® 

10904. 

11289. 

11687. 

12235. 

12666. 

2 

8967. 

9283® 

9718* 

10061® 

10416a 

13904* 

11289. 

11687* 

12235. 

12666. 

3 

8967. 

9283. 

9718® 

10061® 

10416. 

10904* 

11289. 

11687. 

12235. 

12666. 

1 

8322. 

8616® 

9020. 

9338. 

9667. 

10120® 

10477® 

10847. 

11355. 

11756. 

2 

8322. 

8616. 

9020. 

9338. 

9667. 

10120. 

10477. 

10847® 

11355. 

11756* 

3 

8322. 

8616. 

9020® 

9338® 

9667. 

10120® 

10477. 

10847® 

11355. 

11756. 

1 

7726. 

7998® 

8373* 

8669. 

8974® 

9395® 

9726. 

10069* 

10541® 

10913. 

2 

7726* 

7998® 

8373. 

8669* 

8974® 

9395. 

9726. 

10069. 

10541. 

10913. 

3 

7726® 

7998. 

8373. 

8669® 

8974. 

9395. 

9726. 

10069. 

10541. 

10913. 

1 

74505. 

78005. 

82885. 

86778. 

90854* 

96538. 

101073® 

105821. 

112441. 

117723® 

2 

74505. 

78005® 

82885. 

86778. 

90854® 

96538® 

101073. 

105821® 

112441. 

117 723. 

3 

74505® 

78005. 

82885® 

06778* 

90854® 

95533. 

101073* 

135821, 

112441. 

117723. 

1 

30564. 

32000* 

34002. 

35599® 

37271® 

39603. 

41463® 

43411. 

46127. 

48293® 

2 

30564. 

32000. 

34002* 

35599. 

37271. 

39603. 

41463, 

43411. 

46127. 

48293. 

3 

30564. 

32000® 

34002. 

35599. 

37271. 

39603* 

41463® 

43411. 

46127. 

48293. 

1 

24216. 

25354® 

26940. 

28205. 

29530® 

31378® 

32852. 

34395® 

36546* 

38263. 

2 

24216. 

25354. 

26940. 

28205* 

29530® 

31378. 

32852® 

34395. 

36546. 

38263® 

3 

24216® 

25354. 

26940® 

28205. 

29530. 

31373® 

32852. 

34395. 

36546® 

38263, 

1 

20058. 

21000. 

22314® 

23362® 

24459. 

25990* 

27210. 

28489® 

30271. 

31693. 

2 

20058. 

21000* 

22314, 

23362® 

24459. 

25990* 

27210. 

28489. 

30271® 

31693. 

3 

20058. 

21000® 

22314® 

23362. 

24459* 

25993. 

27210® 

28489. 

30271. 

31693. 

<F 

SCll.E-6 

> SC(1. 

E-5) SCU.E-4) 

SC<0.001) 

SC(O.Ol) 

SC(0 

.1) 

SC< 1.0) 

SC(10*0) 


1 

2 

3 


1 * 16525E+D7 
1 * 16525E+07 
i • 16525E+07 


1 ®72128E+07 
1.72128E+07 
1.72128E+07 


2.54265E+07 
2.54265E+07 

2. 54265E+07 


3.75596E+07 
3® 75596E+07 
3, 75596E-9-07 


5® 54823E+07 
5. 54823E+07 
5® 54823E+07 


8 • 1 9575 E+07 
8 » 19575 E+07 
8 » 19575E+07 


la 2 1066E*08 
1* 2 1066E+08 
1. 21066E+08 


1. 78837E+08 
1.78837E*08 
1.78837E+08 


1 

2 

3 


3 » 52454E+04 

3. 52454E+04 
3.52454E+04 


5 .20639E+04 
5.20639E+04 
5 .20639E+04 


7.69079E+Q4 
7 * 69079E+04 
7 • 69079E-9-04 


1 « 13607E+Q5 
1 • 13607E+05 
1 » 13607E+Q5 


1.67B18E+05 
la 67818E+05 
1. 67818E+05 


2.47898E+05 

2.4789BE+05 

2.47898E+05 


3® 66191E+05 
3*66191E+05 
3® 661916+05 


5® 409306+05 
5.40930E+05 
5* 40930E+05 


1 

2 

3 


2.6879OE+03 

2®68790E+03 

2.68790E+03 


3.97052E+03 
3 .97052E+03 
3,970526+03 


5.86519E+03 
5.86519E+03 
5® 86519E+03 


8® 66395E+03 
8.66395E+03 
8® 66395E+03 


1. 27982E+04 
la 27982E+04 
1® 27982E+04 


1 • 89053E+04 
1.89053E+04 
1 « 89053E+04 


2.79266E+04 
2* 79266E+04 
2* 79266E+04 


4* 12527E+04 
4o 12527E+Q4 

4. 1252 7E+04 


1 

2 

3 


8* 75194E+Q2 
8® 75194E+02 
8.75194E+02 


1.29282E+03 
1.29282E+03 
1 .29282E+03 


1.90973E+03 
1.90973E+03 
1 • 90973E+03 


2«82102E+03 
2. 82102E+03 
2.82102E+03 


4® 16717E+03 
4® 16717E+03 
4.16717E+03 


6« 1556SE+03 
6® 15566E+03 
6.15566E+03 


9. 09304E+Q3 
9* 09304E+03 
9® 09304E+03 


1.34321E+04 
1* 3432 1E+04 
i« 3432 1E+04 


1 

2 

3 


3a 16622E+09 
3® 16622E+09 
3o 16622E+09 


5 »90742E+09 
5.90742E+09 
5 ,907426+09 


1® 10218E+10 
1 . 10218E+10 
1.10218E+1Q 


2.05641E+10 

2.05641E+10 

2.05641E+10 


3® 83577E+10 
3« 83677E+10 
3® 83677E+10 


7 • 15849E+10 
7 • 1 5849E+10 
7® 15849E+10 


la 33560E+ 11 
1.33560E+U 
1® 33560E+11 


2. 49191E+ 11 
2 ®49191 E+ 11 
2.49191EH1 


1 

2 

3 


1 • 62858E+04 
la 62858E+04 
1.62858E+04 


3.03854E+04 

3.03854E+04 

3.03854E+04 


5.66919E+04 
5o 6691 9E+04 
5 * 6691 9E+04 


1 • 05774E+05 
1. 05774E+05 
1 « 05774E+Q5 


1® 97348E+05 
1* 973486+05 
1® 97348E+05 


3 * 68204E+05 
3.68204E+05 
3 .68204E+05 


6*869816+05 
6.86981E+05 
6® 86981E + 05 


1*2 81 74 E+ 06 
1 .28174E+06 
1 » 281746+06 


1 2 ® 15354E+0 3 4.01799E+03 7.49660E+03 

2 2® 15354E+03 4.01799E+03 7.49660E+03 

3 2* 153546+03 4.01799E+03 7.49660E+03 


1 « 39869E+04 
1 ®39869E+04 
1.39869E+04 


2.60961E+04 
2* 60961E+04 
2. 60961E+04 


4® 86892E+04 
4*86892E+04 
4.86892E+04 


9.08423E+04 
9* 08423E+04 
9® 08423E+04 


1 ® 69490E+05 
1 .69490E+ 05 
1 • 69490E+05 


1 9.26504E+02 

2 9.26504E +92 

3 9.26504E+02 


lo72863£+03 
1 *72863E+03 
1.72863E+03 


3® 22522E+03 
3.22522E+03 
3* 22522E+03 


6.01748E+03 
6.0174BE+Q3 
6® Q1748E+03 


I® 12272E+04 
1® 12272E+04 
1* 12272E+04 


2® 09472E+04 
2.094726+04 
2 » Q9472E+04 


3.90S25E+04 
3® 90825E+04 
3. 90825E+04 


7 • 29186E+ 04 
7 ® 291 86E+04 
7 » 29186E+ 04 
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EVSIMEAT) 

EVS(CLAD) 

FOOT 

F 

P (GAP ) 

P ( WATER) 

QG/QGMAX 

TCHATER) 

TIME 

0. 

-0. 

-0. 

-0. 

1 • 5000E+01 

1.5000E+31 

-0. 

7 . 3000E+0 1 

1 .OOOOE-O 3 

0. 

-0. 

3.0200E-03 

3.0200E-03 

7.3500E+01 

3.0000E*‘01 

1. OOOOE+OO 

1.7200E+03 

i.OOOOE+OO 


KTl ME= 10 

RADII.... 

FORCES... 

STRAINS.. 

STRESS... 

PARAMS... 

PLASTIC FLOW INDEX * 


DT IME= 

0.97806-01 

LAPSED 

SECOND S= 

5.783 T I ME= 

1.0000 


RA l M 1 — 

0.15130 RB ( M) - 0.31318 

R ( 3 AP ) = 

3 . 54472E-03 RBIC> = 

0.37707 


FCORE= 

0 

FCQN(R)= 0 


FCONI Z ) = 0 



EZ ( M ) = 

0.8615E-02 

EZ{ C )= 0. 5510E-02 

EMAX ( M ) = 0.1565E-03 

EMAX (CI= 

0. 1000E-07 

SMAXIM)= 

8367.7 

SMAXI C ) - 

1546.5 





DV/V = 0.02611 DD/D(C)= 0.00553 DD/D(F-OD>=* 0.00863 DD/0(F-ID)= 0.00864 

111100001100000 


KM KR TEMP RADIUS DEG/DT 


T(AVE) // EV 


1 i 1872.70 0.15130 

-3.23544E-05 0.39998E- 


0. 74825E-03 7438.43 

01 0.39998E-01 0 


4966.26 

0 


1 2 1872.66 0.15432 

-0 • 2 5630E-05 0.39997E- 


0.73475E-03 7405.64 

01 0.39997E-01 0 


5082.51 

0 


1 3 1872.59 0.15634 

-3.24736E-05 0.39998E- 


0.68449E-03 7364.89 

01 0.39998E-01 0 


5041.55 

0 


1 4 1872.22 0.16189 

- 3 • 24140E-05 0.39998E- 


0. 15695E-03 6847.19 

01 0.39998E-01 0 


5423.54 

0 


1- 5 1868.44 0.18348 

-0. 13858E-05 0.39999E 


0. 88742E-05 5891.99 

-01 0.39999E-01 0 


4727.57 

0 


1 6 1851.24 0.20507 

-0 • 491 74E-06 0.40000E 


0. 20562 E-06 4979.48 

-01 0.40000E-01 0 


3274. B4 
0 


1 7 1853.81 0.22676 

-3.11176E-06 0.40000E 


0. 12361E-07 4736.78 

-01 0.40000E-01 0 


1296.60 

0 


1 8 1837.47 0.24835 

-3 • 14901E-07 0.40000E 


0.20037E-05 5489.93 

-01 0. 40000E-01 0 


-1182.18 

0 


1 9 1821.29 0.26993 

-3 • 745Q6E-08 0.40000E 


0 . 2 1662E-03 6980.35 

-01 0.40000E-01 0 


-3751.70 

0 


i 10 1832.38 0.29151 

-3 • 74506E-08 0.40000E- 


0.17172E-03 7664.07 

01 0.40000E-01 0 


-5210.53 

0 


1 11 1783.69 

3 


0.31318 

0 


VH VH(1) VH( 2 ) VHI3) VHI4) VH(5) 


1784.44 

0 0 


1784.39 

0 0 


1784.18 

0 0 


1782.21 

0 0 


1776.99 

0 0 


1768.60 

0 0 


1757.30 

0 0 


1743.26 

0 0 


1726.58 

3 0 


1707.27 

0 0 


0 

0 0 


2 

1 

1743.23 

0.31373 

0 . 16960E- 1 1 

15523.3 

1115.52 

1660.38 




3 

0 

0 

0 

0 

0 

0 

2 

2 

1741.23 

0.31875 

0.21001E-12 

15551.9 

705.898 

1656.66 




3 

0 

0 

0 

0 

0 

0 

2 

3 

1735.41 

0.33384 

0.41637E-15 

15593.6 

112.711 

1651.24 




3 

0 

0 

0 

0 

0 

0 

2 

4 

1729.83 

0.34892 

0* 13366E- 12 

15633.8 

-457.006 

1646.04 




3 

0 

0 

0 

0 

0 

0 

2 

5 

1724.47 

0.36400 

0.1 0404E- 1 1 

15670.1 

-967.711 

1641.36 




3 

0 

0 

0 

0 

0 

0 

2 

6 

1723.00 

0.37707 

0 

0 

0 

0 




3 

0 

0 

0 

0 

0 

0 
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KM KE ET-6 SR SC SZ SG // EG EGR EG£ EGZ EL SL(R-C) SUC-Z) SLIZ-R) 


1 1 23,595 -181,25 -7448,2 -7284,3 7186.4 

0 * 978 IE-04 0.9780E-04 -0.4939E-04 -0.4841E-04 0.9781E-04 0.10D7E+05 -182.45 -9890.3 


1 2 28.595 -378.12 -7465,4 -7419,3 7064.3 

0. 8953E-04 0.8952E-04 -0.44096-04 -0.4543E-04 0,89536-04 9634.4 -20.655 -9613.7 


1 3 28.596 -391.35 -7311.2 -7437.2 6983.7 

0® 8015E-04 0.8010E-04 -0.3798E-04 -0.4212E-04 0.8D15E-04 9171. 1 204.80 -9375.9 


1 4 23,600 -870.58 -7321.1 -8095.5 6870.5 

0.1860E-04 0.1 848E-04 -0.7518E-05 -0.1096E-04 0.1860E-04 6946.2 839.99 -7786.2 


1 5 28.610 -1487*1 -5605.0 -7104.8 5038.1 

0* 9133E— 06 0.8710E-06 -0.2282E-06 -0.6428E-06 0.9133E-06 4138.9 1507.7 -5646.6 


1 6 28.627 -1769,7 -3132.7 -4931.9 2747.2 

0.3112E-07 0.1721E-07 0.1990E-08 -0.1920E-07 0.3112E-07 1363.3 1799,6 -3162.9 


1 7 28,649 -1754.4 -174.52 -1964.8 1694.9 

0.1124E-07 -0.5272E-09 0.1238E-08 -0.7110E-09 0.U24E-07 -1579.9 1790.3 -210.37 


1 8 23.677 -1438.9 3252.3 1736.6 4146.6 

0. 210 IE-06 -0.1899E-06 0.1496E-06 0.4026E-07 0.2101E-06 -4697.7 1517.8 3179.9 


1 9 2B.711 -951.46 6550.3 5667.5 7101.6 

0* 2493E-04 -0.2475E-04 0.1488E-04 0.9867E-05 0.2493E-04 -8260.2 978.70 7281.5 


1 10 28.749 -362.71 8004.9 8005.0 8367.7 

0.1565E-03 — 0 • 1 564E-03 0.8198E-04 0.7440E-04 O.I565E-03 -0.1294E+05 145.29 0.1279E+05 


2 1 23.286 -89.378 -1553. V -1706.6 1546.5 

0. 1000E-07 0.1668E-12 -0.7040E-13 -0.9644E-13 0.1000E-07 1464.5 152.70 -1617.2 


2 2 23,292 -115.89 -910.42 -1093.5 900,14 

0. 10G0E-Q7 0.2022E-13 -0.6990E-14 -0.1323E-13 0.1000E-07 794.53 183.07 -977.61 


2 3 23.301 -130,83 -4,7111 -202,93 173.77 

0. 1000E-07 0 0-0 0. 1000E— 07 -126.12 198.22 -72.093 


2 4 23.310 -105.79 826.78 651.40 858.42 

0. 1000E-07 — 0.1287E-13 0.8436E-14 0.4429E-14 0.1000E-07 -932.56 175,37 757.19 


2 5 23.318 -57.658 1545.5 1418.2 1543.4 

0» 1000E-07 -0*102 5E-12 0.5703E-13 0.4546E-13 0.1000E-07 -1603.2 127.30 1475.9 


KM 

KE 

KG 

DFT 

PGAS 

PA 

SG(AVE) 

SP 

SC 

RA 

V A 

VPC 

NH( T ) 

MOLES 

1 

1 

1 

1 » OOE+OO 

6.73E+01 

-8.01E+03 

6.05E+03 

8.65E+03 

1.13E+08 

7.07E-07 

4.00E-02 

-1.82E-01 

2.70E+16 

6.47E-08 

l 

2 

1 

l.OOE+OO 

6.73E+01 

—8.0 1E+Q3 

6.10E+03 

8. 65E+03 

1.08E+08 

7.07E-D7 

4.00E-02 

-1.90E-01 

2.70E+16 

6. 47E-08 

1 

3 

1 

1 , OOE+OO 

6.73E+01 

-8.01E+03 

6.08E+03 

8.65E+03 

1 * 10E+08 

7* 07E-07 

4.00E-02 

-1.87E-01 

2.70E+16 

6.47E-08 

1 

4 

1 

l.OOE+OO 

6.73E+01 

-8.01E+03 

6.26E+03 

8.65E+03 

9.89E+07 

7.07E-07 

4.00E-02 

-1*7 9E-01 

2. 70S + 16 

6.47E-08 

1 

5 

1 

l.OOE+OO 

6.71E+01 

-8.01E+03 

5.93E+03 

8, 66E+03 

1.43E+Q8 

7.07E-07 

4.006-02 

-1.50E-01 

2* 70E + 16 

6. 47E-08 

1 

6 

1 

1 * OOE+OO 

6.69E+01 

-8.01E+03 

5.26E+Q3 

8.66E+03 

2.92E+08 

7.07E-07 

4.00E-02 

-1.12E-01 

2,706+16 

6.47E-08 

1 

7 

1 

l.OOE+OO 

6.65E+01 

-8.01E+03 

4.346+03 

8.67E+03 

9.10E+08 

7.07E-07 

4.006-02 

-7.64E-02 

2.70E+16 

6.47E-08 

1 

8 

1 

l.OOE+OO 

6.61E+01 

-8.01E+03 

3.18E+03 

8.68E+03 

5 * 266+09 

7.07E-07 

4.00E-02 

-4.82E-02 

2.70E+16 

6.47E-08 

1 

9 

i 

l.OOE+OO 

6.56E+01 

— 8.01E+03 

1.98E+03 

8 .696+03 

7.06E+10 

7.07E-07 

4.00E-02 

-2* 7 IE-02 

2.70E+16 

6.47E-08 

1 

10 

X 

l.OOE+OO 

6* 50E+01 

-8.01E+03 

1.30E+03 

8.71E+03 

7.64E+11 

7 . 07E— 07 

4.00E-02 

-2.83E-02 

2.70E+16 

6.47E-08 


0 KRDT 

KCTAC 

JERKS 

RQPLUG 

0 

10 

0 

l.OOOOE+OO 
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EVS( MEAT ) EVStCLAD) 
4. 7000E-04 -0* 

1. 1800E-02 -0. 


FDQT F 

3 . 0200E-03 3.0200E 

3 . 0200E— 03 7.5500E 


P l G AP ) P 

00 7.9400E+01 3 

01 2.2100E+02 3 


WATER) QG/QGMAX 
OOOOE+Ol 1 „ OOOOE+OO 
OOOOE+Q l l.OOOOE+QO 


T (WATER ) TIME 
1 . 7200E+0 3 1.0000E+03 

1.7200E+03 2. 50006+04 


KTIME* 


53 

DT IME* 

316.46 

LAPSED 

SECONDS* 41.317 

TIME* 18906.9277 

RADII... 

• 


RA(M)= 

0.15156 

RB l M) = 0.31373 

R { 3 AP ) = 0 

RBIC )= 0.37708 

FORCES.. 

• 


FCORE* 

0 

FCON ( R ) = 0 

FCONt Z ) = 

0 

STRAINS. 

- 


EZ (M ) = 

0.1041E-01 

EZ ( C ) = 0. 5514E-02 EMAX(M)= 

0.4666E-03 EMAX(C)= 

STRESS.. 

* 


SMAX(M) 

* 652.13 

SMAX ( C ) = 

1963.2 


PARAMS.. 

PLASTIC 

FLOW 

INDEX = 

DV/V = 
0 0 0 0 

0.03154 
0 0 0 0 0 0 0 

DD/D (C ) = 0,00555 
0 0 0 0 

OD/D ( F-OD ) = 0.01041 

D D/D ( F- I D ) = 0.01042 


EVS ( MEAT ) EVSICLAD) FDOT F P(GAP) P ( WATER ) QG/QGMAX TtWATER) TIME 

2® 3500E-02 -0„ 3.0„200E-03 1.5100E+D2 3.6800E+02 3.0000E*Ql l.OOOOE+OO 1.7200E+03 5.0000E+04 

2® 3500E-02 -0. -0. 1.5100E+02 7.5000E+01 1.5000E+D1 -0. 7.0000E+01 5.0D01E+04 

KT I ME= 833 DT IME= 0.4501E-01 LAPSED SECONDS* 618.100 TIME* 50001.0000 

RADII..®. RA (M ) = 0.14000 R8(M)= 0.3L172 R(5AP)= 0.67439E-03 RB(C)= 0.37532 

FORCES... FCOR E = 0 FCON(R)= 0 FCON(Z)= 0 

STRAINS.. EZ(M)= 0.2489E-02 EZ{C)= 0.1803E-03 EMAX(M)= 0.1125 EMAX(C)= 0.1588E-02 

STRESS... SMAX ( M ) = 0.1405E+05 SMAX(C)= 1517.7 

PARAMS. .. DV/V = 0.05217 DD/DIC)* 0.00084 DD/D(F-OD)= 0.00394 DD/ D ( F- ID > =-0 .06669 

PLASTIC FLOW INDEX = 010000000100000 


KM KR TEMP RADIUS DEG/DT SP PB T(AVE) // EV VH VH(1) VH(2) VH(31 VH(4) VH(5) 


1 1 70.0002 0.14000 0.50514E-22 12147.7 -8465.04 111-940 

0.60082E-01 0.76582E— 01 0.76582E-01 0 0 0- 0 


1 2 70.0002 0.14339 Q.12879E-21 12125.5 -8980.41 111.938 

0.59721E-01 0.76221E-01 0.76221E-01 0000 


1 3 70.0002 0.14565 0.23643E-22 12093.6 -8072.34 111.929 

0 • 59538E-01 0.76038E-01 0.76038E-01 0000 


1 4 70.0002 0.15182 0.72291E-23 11985.9 -7589.60 111.847 

0.57620E-01 0.74 120E-0 l 0.74120E-01 0000 


1 5 70.0002 0.17546 0.49574E-24 11834.2 -6074.03 111.635 

0. 5 5954E-01 0.72454E-01 0.72454E-01 0 0 0 0 


1 6 70.0002 0.19869 0.11004E-25 11699.6 -4012.85 111.306 

0.54021E-01 0.70521E-01 0.70521E-01 0000 


1 7 70.0002 0,22172 0.90928E-26 11577.8 -1448.21 110.873 

0.51950E-01 0 . 68450E-01 0.68450E-01 0000 


1 8 70.0002 0.24443 0.56144E-24 11466.6 1570,14 110.346 

3.49921E-01 0 .6642 IE— 01 0.66421E-01 0 0 0 0 


1 9 70.0001 0.26696 0.10178E-22 11362.3 5022.48 109.730 

0.47847E— 01 0.64347E-01 0.64347E-01 0000 


1 10 70.0001 0.28935 0.54126E-20 11263.2 8901.66 109,029 

0.4580 5E-01 0 « 62 30 5E-01 0.62305E-01 0 0 0 0 


1 11 70,0001 0.31172 0000 

0 0 0 0 0 0 0 
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2 

1 

70,3000 

0.31240 

0. 14973E-28 

73992,5 

618.598 

107.788 




3 

0 

0 

0 

0 

0 

0 

2 

2 

70,3000 

0.31738 

0. 11157E-29 

73618.2 

366.517 

107.673 




3 

0 

0 

0 

0 

0 

0 

2 

3 

70,3000 

0.33236 

0 

73065.0 

-6.94080 

107.505 




0 

0 

0 

0 

0 

0 

0 

2 

4 

70,3000 

0.34735 

0.663976-29 

72526.7 

“372.738 

107.343 




3 

0 

0 

0 

0 

0 

0 

2 

5 

70# 3000 

0.36233 

0.38927E-28 

72033.9 

“711.012 

107.199 




3 

0 

0 

0 

0 

0 

0 

2 

6 

70.3000 

0.37532 

0 

0 

0 

0 




3 

0 

0 

0 

0 

0 

0 


KM 

KE 

ET-6 SR 

sc SZ 

SG // 

EG EGR 

EGC EGZ 

EL SL{ R-C 1 

SL(C“Z1 

SHZ“RI 

1 

1 

33.938 

81.373 

0.1299E+05 

0. 1235E+05 

0* 1260E+05 






0.1125 

0.1043 

“0.8687E-01 

“0.1746E-01 

0.1118 

0.3930E+07 

“0. 1431E+07 

“0.249 9E+07 

1 

2 

33.938 

339.67 

0. 1391E+05 

0. 1271E+05 

0.1302E+05 






0. 1087 

0.1008 

“0.8341E-01 

-0. 1735E-01 

0.1377 

0.378 5E+07 

-0. 1361E+07 

-0. 2423E+0 7 

1 

3 

33.938 

677.91 

0.116BE+05 

0* 1188E+05 

0.1111E+05 






0.1027 

0. 9588E-01 

-0.7860E-01 

“0. 1728E-01 

0.1022 

0.3588E+07 

-0.1265E+07 

-0. 2323E+07 

1 

4 

33.938 

1565.4 

0 . 1009E+05 

0. 11 14E+05 

9092.9 






0.8558E“01 

0. 809 1E“01 

“0.6426E-01 

-0.1666E-01 

0.8545E-01 

0.2986E+07 

-0.9828E+06 

—0. 2003E+07 

1 

5 

33.938 

2429.0 

6801.6 

9009.7 

5800.9 






0. 6629E-01 

0 * 6367E— 01 

“0.4758E“01 

-0.1609E-01 

0 . 6621E-01 

0.2290E+07 

“0.6518E+06 

-0.1638E+07 

1 

6 

33.938 

2718.7 

3307.5 

6024.4 

3054.1 






0. 5270E-01 

0* 5130E-01 

-0.3589E-01 

-0. 1541E-01 

0.5265E-01 

0. 1798E+07 

-0.4252E+06 

-0. 1373E+07 

1 

7 

33.939 

2591.1 

-498.95 

2256.8 

2937.2 






0.4272E-01 

0.4204E— 01 

-0.2737E-01 

-0.1467E-01 

Q.4267E—01 

0. 1435E+07 

-0.26486+06 

-0.1 170E+07 

1 

8 

33.939 

2135.4 

-4637.9 

-2212.7 

5944.1 






0. 3532E-01 

0.3491E“01 

“0.2100E-01 

-0. 1392E-01 

0.3515E-01 

0.1160E+07 

“0. 1484E+06 

-0.1012E+07 

1 

9 

33.940 

1407.9 

-9124.4 

-7366.0 

9772.5 






0. 2974E-01 

0.2920E“01 

-0. 1606E-01 

“0* 1314E“01 

0.2925E— 01 

0. 9441E+06 

-0.6207E+05 

-0.8820E+06 

1 

10 

33.940 

442.20 

“0. 1399E*05 

-0.1318E+05 

0. 1405E+05 






0* 2539E-01 

0.2452E— 01 

“0. 1217E“01 

“0. 1235E-01 

0.2452E-01 

0. 7712E+06 

2898.7 

-0.7741E+06 


2 

1 

25.874 -84.863 

0® 1588E-02 -0* 1452E— 02 

-1451.4 

0,12836-02 

-321,39 

0* 1687E-03 

1265.0 

0® 1588E— 02 

-0.4580E+05 

0.1809E+05 

0,277 1E+05 

2 

2 

25.874 -113.52 

0.1470E-02 -0* 1 350E—02 

-801.40 

0® 1177E-02 

-185.73 

0. 1728E-03 

654.77 

0. 1469E— 02 

-0.4291E+O5 

0.1671E+05 

0.262 0E+05 

2 

3 

25.874 -124,67 

0.1310E-02 -0 ® 1212E—02 

112,20 

0® 1035E-02 

33.311 

0.1778E-03 

208.91 

0. 1310E— 02 

-0.3900E+05 

Os 1486E+05 

0, 241 4E+05 

2 

4 

25.874 -97,070 

0* 1 170E-02 “0. 1092E—02 

947.44 

0,91016-03 

268.96 

0, 1815E-03 

917.96 

0, 11 70E-02 

-0.3557E+05 

0, 1325E+05 

0* 2232E+05 

2 

5 

25.875 -44.687 

0* 10556-02 -0.9908E-03 

1671 * 9 

0.8072E-03 

507,98 

0. 1836E-03 

1517.7 

0.1054E-02 

-0® 3273E+05 

0.1192E4-05 

0.2081E+05 
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KM 

KE 

KG 

3FT 

PGAS 

PA 

SG(AVE) 

SP 

SC 

RA 

VA 

UPC 

NHIT) 

MOLES 

1 

1 

1 

l.OOE+OO 

1.26E+03 

-5.24E+03 

1 » 86E +03 

9.93E+03 

2o30E+30 

8.78E-07 

7« 66E-02 

4*0 0E-02 

2® 70 E + 16 

2® 06E-05 

1 

2 

1 

l.OOE+OO 

1 » 27E + 03 

-5.25E+03 

2.15E+03 

9.93E+03 

1 * 13E + 30 

8.77E-07 

7.62E-02 

4*6 IE— 02 

2.70E+16 

2® 06E-05 

1 

3 

1 

l.OOE+OO 

1.27E+03 

-5®25E+03 

1 o 62E+03 

9.93E+03 

4«47E+30 

8.76E-07 

7® 60E-02 

3.47E-02 

2® 70E+1 6 

2.06E-05 

1 

4 

1 

l.OOE+OO 

1 c29E+03 

-5.28E+03 

1. 31 E +03 

9.93E+03 

1.27E+31 

8.69E-07 

7.41E-02 

2.75E-02 

2.70E+16 

2® 06E-05 

1 

5 

1 

l.OOE+OO 

1.31E+03 

-5.31E+03 

4.29E+02 

9.93E+03 

3.07E+33 

8.62E-07 

7.25E-02 

8* 73E-03 

2® 70E+16 

2.06E-05 

1 

6 

1 

l.OOE+OO 

1.33E+03 

-5.35E+03 

7.47E+02 

9 * 93E+03 

2,04E+32 

8.54E-07 

7.05E-02 

-1.54E-02 

2* 70E + 16 

2® 06E-05 

1 

7 

1 

l.OOE+OO 

1.37E+03 

— 5* 38E+03 

2.18E+03 

9.94E+03 

1 * 09E+30 

8» 46E-07 

6# 85E— 02 

-4.36E-02 

2* 70 E+l 6 

2.06E-05 

1 

8 

1 

l.OOE+OO 

1.40E+03 

-5.41E+03 

3« 83E + 03 

9.94E+03 

7.06E+28 

8. 37 E— 07 

6® 64E-Q2 

-7.54E-02 

2® 70E + 16 

2* 06E-05 

1 

9 

1 

l.OOE+OO 

1.45E+03 

— 5«44E+03 

5* 67E+03 

9.94E+03 

1.Q5E+28 

8„ 29E— 07 

6.43E-32 

-1.30E-01 

2.70E+16 

2* 06E-05 

1 

10 

1 

1 . OOE+OO 

1.51E+03 

-5.46E+03 

7.70E+03 

9.94E+03 

2.41E+27 

8® 20E-07 

6.23E-02 

-2* 56E-01 

2.70E+16 

2.06E-05 


0 

KRDT 

KCTAC 

JERKS 

RQPLUG 


32 

1131 

0 

1 * QOOOE-Q6 


EVS ( MEAT ) EVS ( CLAD) FOOT F P (GAP ) P ( WATER) QG/QGMAX T (WATER ) TIME 

2.350QE-02 -0, -0® U5100E + 02 7.5000E+01 1.5000E+01 -0* 7«OOOOE+Ol 5*0001E+04 

BLANK CARD END OF PROBLEM 

KRDT KCTAC KTIME TINY DELTA JERKS RDT STEPS/SEC., MINUTES 

32 1131 833 1 « OOOOE-07 4.5010E-02 0 l.OOOOE+OO 1.3442E+00 1.0328E+01 


END OF OUTPUT 

»01* EXIT IN SKEN 
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APPENDIX C 


Subroutine 

SWELL 

STIM 

STEM 

SEND 

SGAS 

SKEN 

SRAW 

LATON 

LATON1 

FINIS 

LUST 

SEAS 


CTAC 

SEGT 


MODIFIED CYGRO-2 SUBROUTINES AND MAJOR FUNCTIONS 

Function 

main program: controls all input data card reading 

main calculational program: controls all calculations, selects time incre- 
ments, initializes, and prints output 

calculates heat generation and temperature distributions 

sets up interpolation of time cards and controls interpolation of properties 
table 

gas bubble calculations, solid swelling 

reads time cards, input transition to other subroutines 

input transition between reading subroutine and storing subroutine 

reads all input data 

former part of subroutine LATON 

writes 'end of output* 

interpolation procedure for tables 

finishes calculating stress and strain increments and adds to value at start 
of time step, interpolates material properties as a function of tempera- 
ture, prints some output, performs matrix manipulations, and serves a 
transition for gas bubble calculations 

calculates contact loads between rings, performs matrix manipulations 

calculates creep rates and plastic flow rates, interpolates properties as a 
function of flux or flux rate 
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TABLE I. - PROPERTIES OF T-lll 


[Poisson's ratio, 0. 35. ] 


Temperature 

Modulus of elasticity 

Coefficient of expansion 

Conductivity 

°F 

K 

psi 

N/cm 2 

in. / (in. )(°F) 

cm/ (cm)(K) 

Btu/ (hr) (in. )(°F) 

W/(cm)(K) 

32 

273 

26.0X10 6 

18. 2xl(j® 

3-OOxlO -6 

5.40X10 -6 

2.00 

0.415 

1832 

1273 

23.0 

16.1 

3.65 

6.57 

2.81 

.583 

2400 

1589 

22.0 

15.4 

3.92 

7.06 

3.00 

.623 

2732 

1773 

21.4 

15.0 

4.08 

7.34 

3.05 

.633 


TABLE II. - PROPERTIES OF UN 
[Poisson's ratio, 0. 28. ] 


Temperature 

Modulus of elasticity 

Coefficient of expansion 

Conductivity 

°F 

K 

psi 

N/cm 2 

in. / (in. )(°F) 

cm/(cm)(K) 

Btu/ (hr) (in. )(°F) 

W/(cm}(K) 

32 

273 

31.0X10 6 

21.7X10 6 

4. 67X10" 6 

8.40X10 -6 

0.770 

0. 160 

932 

773 

30.3 

21.2 

4.67 

8.40 

. 982 

.204 

1832 

1273 

28.5 

19.9 

5.06 

9. 11 

1.160 

.241 

2732 

1773 

22.8 

15.9 

5.34 

9.61 

1.248 

.259 


40 





TABLE HI. - STRESS- STRAIN AND CREEP 


EXPONENTS OF T- 111 


Temperature 

Stress-strain exponents 

Creep exponents 

°F 

K 

P 

Q 

C 

G 

32 

273 

-31.270 

15.10 

-30.000 

3.692 

1832 

1273 

-25.427 



-10.474 



2400 

1589 

-23. 900 



-7.230 



2732 

1773 

-22.665 

1 


-5.878 




TABLE IV. - STRESS-STRAIN AND CREEP 
EXPONENTS FOR UN 


Temperature 

Stress-strain exponents 

Creep exponents 

°F 

K 

P 

Q 

C 

G 

32 

273 

-22.699 

20.0 

-30.000 

5.902 

932 

773 

-22.053 



-15.131 



1832 

1273 

-21.405 



-8. 534 



1732 

1773 

-20.759 



-5.658 







TABLE V. - INPUT DATA CARDS 


Card 

Comments 

Comment cards 

These are the first cards to appear after the $DATA card. The modified CYGRO-2 
allows any number of comments on as many cards as are necessary to completely 
describe the problem. These words can fill columns 3 to 72. The last comment card 
must contain a numerical number greater than 0 in columns 1 to 2. The next card must 
be the 900 card. 

900 

The first variable on the 900 card is the plastic flow rate exponent. Allowing the expo- 
nent to approach plus infinity would make flow instantaneous. However, this is not 
practicable on the computer. On the IBM 7094 an exponent is limited in magnitude and 
it was found that a value of 30 to 35 was satisfactory. This value allowed sufficiently 
rapid relaxation of stresses, whereas larger exponents occasionally created run 
stoppage . 


The run limit specified in this version of CYGRO-2 was and is not utilized. It is simply 
erased when the next card is read. 


In all tabulated input data, only the range between maximum and minimum specified 
values is used. If at any time the tabulated variable has an incompatible range as com- 
pared with a calculated variable, the end points are used; no extrapolation beyond the 
table ever occurs. 

90K 

K = 1, 2, 3 

The gas leakage function is difficult to evaluate for problems where cycling takes place. 
The function on page 19 reference 1 must be integrated over the life of the fuel pin to 
determine the correct coefficient or function if a specified amount of gas release is 
required. 

904 

The film coefficient between fuel and clad can create problems if the gap coefficient and 
no- gap coefficient are greatly different. Since the temperature level in the fuel changes 
suddenly upon contact (it gets cooler), the fuel would thermally contract away from the 
clad. An oscillation can be set up where the program stays in the "twilight zone” for 
a period of time. 

905 

The maximum allowable stress and strain change (AS and Ae, respectively) ^are not 
computer-time restrictive based solely on magnitude. There is some optimum AS 
and Ae for each problem. For example, a problem ran with Ac = 0. 001 in. /in. 

(0. 001 cm/cm) and AS = 1000, 2000, and 10 000 psi (689, 1379, and 6895 N/cm^) ran 
the quickest with AS = 2000 psi (1379 N/cm ). Too large an increment can set up in- 
stabilities which take longer to straighten out. Maximum allowable temperature in the 
plug must be approached with caution. The CYGRO-2 program calculates a tempera- 
ture distribution; if the maximum is exceeded, the volumetric heat generation rate 
specified on card 1000 is proportioned down such that the centerline or plug tempera- 
ture is exactly the specified maximum. There may be occasions when this is not 
desired. 
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TABLE V. - CONTINUED. INPUT DATA CARDS 


Card 

Comments 

906 

The calculation option, Sg m * n , on the form of partial derivatives represents a min- 

imum stress at which this option is utilized. As a stress, it is not critical that it 
equal 1. Only a small number is required for the same result. 

Fuel-plug contact requirement is necessary for the case of solid fuel pins. It is this 
variable which allows the user to specify a solid or hollow fuel pin in addition to radii 
specification. 

Input- output option, Mg ag , determines whether to expect only one 1002 card or a card 
for each individual ring. Some of the data specified on more than one 1002 card are not 
stored; these are discussed under the 1002 card comments. The E^^g option for 
solid fuel swelling must be used with the pressure factor R vp on card I0J0. This 
allows the user to lump all sources of swelling into one modified solid swelling if de- 
sired. This approach to modeling fuel swelling may be desirable if data on swelling as 
a function of pressure and temperature are available. 

907 

Two new variables have been added to this card. These variables are counters dealing 
with the printout of the calculations and appear in the seventh and eighth data fields, 
columns 57 to 72. They are print frequencies. Each is given a number, m and n, 

ii. 

respectively. Every m time step a printout of group 4 (ref. 1) data occurs, and 
every n th step a printout of groups 4 to 7 occurs. This feature easily allows more of 
the calculations to be printed for cases where long periods of time occur between his- 
tory cards. 

1000 
I = 1, 2 

i 

! 

Two factors or multipliers are used to modify the material mechanical properties 
around the gas bubbles. One is a creep stress factor and the other is a plastic stress 
factor; they appear in the third and fourth fields of the card. The plastic stress factor 
appears as a multiplier on the 0. 2 percent yield stress. Since, for bubble equilibrium 
calculations^ strain hardening is not considered, this provides an alternate method of 
inputing some average strain hardening. The creep stress factor is a multiplier also 
utilized in determining the stress level around the bubble. For bubbles calculations 
the stress is never allowed to exceed the yield; and if yield is reached, the extent to 
which plastic flow extends into the material region around the bubbles must be deter- 
mined. A plastic stress factor greater than 1 increases the yield point and, conse- 
quently decreases the plastic zone (although the average stress may be higher). A 
creep stress factor greater than 1 increases the plastic zone. It also increases the 
apparent stress level in the creep zone. These factors can be arbitrarily used to in- 
crease or decrease fuel swelling. However, for UN fuel, use of factors equal to 1 has 
resulted in satisfactory agreement between experimental and calculated fuel swelling. 

Full- power heat generation rate must be the reference power level if the fuel pin is 
cycled. 

1002 
I = 1, 2 

These cards describe the gas bubble conditions in the fuel and clad. If all the rings are 
the same, there need to be only as many cards as there are bubble types. Some of these 
conditions may be varied from ring to ring. If this is done, there should be a card for 
each ring. 
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TABLE V. - CONCLUDED. INPUT DATA CARDS 


Card 

Comments 

1002 (continued) 

The fission coefficient was subscripted in the modified CYGRO-2 such that it may also 
vary from ring to ring. This will allow, in effect, a variable gas generation rate in the 
radial direction. Total solid volume associated with each bubble type is not ring depend- 
ent within the program. This is a fractional value, and in most cases should total up to 
1.0. However, it is not mathematically necessary. 

1004 
1 = 1, 2 

Gas pressures used in this tabulation must be specified and they must be in three dis- 
tinctly different values. Even if the pressure factor is not utilized, the three speci- 
fied pressures cannot all be zero or the same. 

1005 
1 = 1, 2 

Any number of node points, from 2 to 4, may be specified for a constant or variable 
heat generation, as a function of the radius index. However, the complete range of the 
index must be covered in nondecreasing order. 

IOJK 
I = 1, 2 
J = l, 2, 3, 4 
K = 1, 2, 3 

It is characteristic of most materials that the creep coefficient decreases (in absolute 
value) nonlinear ly with increasing temperature while the exponent remains constant or 
nearly constant. Interpolation between two temperatures in the table always results in 
a creep exponent that is too large. The result is an apparently stronger material. It 
can appear stronger by a few orders in magnitude in creep rate. This situation, it has 
been found, is best alleviated by setting up the table around an operating temperature 
or temperature range. In this way the range over which interpolation occurs is much 
smaller, and the resulting error is also much smaller. Caution must be used, in any 
case, to prevent the operating temperature from falling midpoint between two widely 
separated temperatures in the table. An example of spacing the temperatures in the 
table in a more useful manner is given in the sample problem for the cladding. 


Another interpolation problem was encountered in the temperature-dependent plastic 
stress-strain exponents. It at first seemed appropriate to solve for the exponents at 
each temperature in the table by forcing the function through two points - the 0. 2 per- 
cent yield point and another point at much higher plastic strain which would be repre- 
sentative of maximum strains expected for the problem. However, the resulting coef- 
ficients and exponents usually behave rather badly; that is, they do not plot as a smooth 
function with temperature. The resulting interpolation can generate a pair of coeffi- 
cients and exponents which do not correctly model the material behavior, but can result 
in a 0.2 percent yield, for example, which is outside the range of intended interpola- 
tion. Sometimes this can cause an instability in the stresses and strains which results 
in termination of the computer problem. It has been found that the best results are 
obtained when the following procedure is used: First, obtain the coefficients at the ex- 
pected operating temperature; then, keeping the stress exponent constant or nearly so, 
vary the coefficient such that the function will pass through the 0. 2 percent yield point 
at each temperature level. Keeping the stress exponent constant and only varying the 
coefficient results in a curve showing more strain hardening at lower temperatures 
than at higher temperatures, which is usually the way most materials behave. If, how- 
ever, it is desired to vary both coefficients and exponent to properly model the mater- 
ial, caution should be exercised to avoid instabilities in interpolation. 


The first data field utilized on this card is the second field. The first field is left blank. 
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Temperature 

°F(K) 



Strain, in. /in. (cm/cm) 

Figure 5. - Stress-strain diagram forT-111 alloy. 





Figure 6. - Stress-strain diagram for uranium nitride. 


NASA- Langley, 1971 
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